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ABSTRACT 

We describe a numerical implementation of star formation in disk galaxies, in which 
the conversion of cooling gas to stars in the multiphase interstellar medium is gov- 
erned by the rate at which molecular clouds are formed and destroyed. In the model, 
clouds form from thermally unstable ambient gas and get destroyed by feedback from 
massive stars and thermal conduction. Feedback in the ambient phase cycles gas into 
a hot galactic fountain or wind. We model the ambient gas hydrodynamically using 
smoothed particle hydrodynamics (SPH). However, we cannot resolve the Jeans mass 
in the cold and dense molecular gas and, therefore, represent the cloud phase with 
ballistic particles that coagulate when colliding. We show that this naturally produces 
a multiphase medium with cold clouds, a warm disk, hot supernova bubbles and a hot, 
tenuous halo. Our implementation of this model is based on the Gadget N-Body code. 
We illustrate the model by evolving an isolated Milky Way-like galaxy and study the 
properties of a disk formed in a rotating spherical collapse. Many observed properties 
of disk galaxies are reproduced well, including the molecular cloud mass spectrum, the 
molecular fraction as a function of radius, the Schmidt law, the stellar density profile 
and the appearance of a galactic fountain. 
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1 INTRODUCTION 

Galaxies form when gas cools radiatively inside a dark mat- 
ter (DM) halo. The haloes in turn form from the non-linear 
collapse of initially small density perturbations that were 
imprinted by quantum fluctuations before inflation. When 
the virial temperat ure of the halo is high e nough the gas 
can co ol radiatively dRees fc Ostriked (|l977l h lwhite fc Reesl 
(1978)) and may become self-gravitating. Some fraction of 
the cooling gas can form stars, which then affect the baryon 
distribution and star formation rate through feedback; for 
examp le from supernova (SN) explosions (e.g. iDekel fc Silkl 

Simulating the growth of dark matter haloes from ini- 
tially small cosmological density perturbations has become 
routine to the extent that even the complex n on-linear stage 
can be predicted with relative confidence fe.g. lSpringel et al 
(2005)). However, following the behaviour of the baryons un- 
til stars form is much more demanding for two main reasons. 
Firstly, the physical processes that lead to star formation are 
still relatively poorly understood, even in the Milky Way 
(MW). Stars are thought to form in molecular clouds in a 
complex interstellar medium (ISM) in which magnetic fields 
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Safi er. McKee fc StahleU (119971 ) ) , osmic ray s, turbulence 
Krumholz fc McKeel (|2005l )h relativistic jets l|Klamer et all 
molecules, dust, and radiative transfer may all play 
some role. Secondly, the physical scales on which star forma- 
tion takes place is vastly different from those of cosmological 
interest. Therefore, simulations on a sufficiently large scale 
to sample cosmological structures cannot presently also re- 
solve physics on the scales relevant for star formation. For 
these two reasons simulations invariably include 'sub-grid' 
physics that model the complexity of star formation in the 
interstellar medium using simple rules. 

In the first generation of hydrodynam ic mod- 



els for galaxy f or mation dNavarro & White 
Steinmetz fc Mullerl ll 1994) : ICen fc Ostrikerl 



1993) 



1992) 



Katz. Weinberg fc HernquistJ (11996I D. the gaseous com- 



ponent of galaxies is modelled as a single fluid. Because 
the spatial and mass resolution in these simulations is 
insufficient to resolve star formation they rely upon a 
simple star formation prescription. In its simplest form this 
consists of devising criteria by which simulated gas can 
be flagged as eligible to form sta rs and then conv erting 
the gas to stars by hand (see e.g. iKav et all (|2002r i for a 
comparison of such criteria). Feedback from core collapse 
(type II) SNe, modelled as an extra source of thermal or 
kinetic energy, was found to have little effect in the early 
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models. This is because the gas in the surroundings of star 
formation sites is at high density and so is very efficient at 
just radiating away the added energy. As a consequence, too 
much gas cooled in dense knots, producin g galaxies much 
more c oncentrated than those obser ved dNavarro fe Benj 
ljl99ll ): IWeil. Eke fc Efstathioul <|l998h >. 

The ISM in observed galaxies is much more complex 
than the single phase medium present in these early sim- 
ulations. Although observed galaxies may have comparable 
mean gas density to the simulated ones, the real ISM consists 
of very dense cold clouds with small volume filling factor, 
embedded in a much more tenuous low-density, hot medium 
with a 'warm' phase at the interface. SN explosions in the 
tenuous medium have a much greater impact on the galaxy 
because this medium cools much less efficiently than gas at 
the mean density. Throughout this work we will use the la- 
bel 'cold' to describe the molecular clouds in the ISM; and 
the labels 'warm' and 'hot' to describe the properties of the 
ambient gas phase at temperatures of approximately I0 4 K 
and 10 6 K respectively. 

In response to these problems with the simplest star 
formation and feedback criteria several authors have intro- 
duced 'multiphase' models for star formation in which the 
ISM is treated as a number of distinct phases. These schemes 
take variou s forms including m o difica ti on of the simulatio n 
algorithm ilRitchie fc Thomas! l|200ll) ; ICroft et all (|200Cj ): 
IScannapieco et al l 20061 )'). treating the multiphase medium 
implicitly by formulating differential e quations tha t mode l 
the interactions between t h e phases (|Yepes et al.l l|l997h ; 
ISpringel fc Hernquistl (|2003l ) ; lOkamoto et all 1120051 )) , explic- 
itly decouplin g supe rnova heated gas from its surroundings 
l|Stinson et all {2006)), or by decoupling the cold molecu- 
lar phase from the hot phase b y one of a variety o f meth - 
ods, including 'sticky pa r ticles' l|Semelin fc Combei (j2002); 
lHarfst. Theis fc Henslerl (|2006j))_ or removing 'cold ' part i- 
cles from the SPH calcula t ion ( Hultman fc Pharasvnl l|l999l) ; 
iPearce et all (| 19991 . 1200 ll ); iMarri fc White! (|2003l )) The de- 
coupling of the hot and cold ISM phases allows thermal heat- 
ing from SN feedback to become more efficient (due to the 
much lower density of the hot phase) and also allows one to 
follow the properties of the cold molecular phase of the ISM. 

This paper describes an attempt to mimic the mul- 
tiphase medium in a star forming galaxy. Since stars are 
known to form in molecular clouds our method focuses on 
simple rules for cloud formation. The star formation recipe 
then simply converts the most massive of these clouds, 'Gi- 
ant Molecular Clouds' (GMCs), into stars with an imposed 
efficiency taken from MW observations. Once a GMC forms 
stars stellar feedback destroys the cloud. We do not attempt 
to model the clouds themselves using hydrodynamics, be- 
cause our current galaxy formation simulations lack by some 
margin the dynamic range to resolve the Jeans mass (Mj) 
in the cold, dense gas that makes up the clouds. We can 
demonstrate this lack of re solution by conside ring theoreti- 
cal models of the ISM (e.g. lWolfire et all (|l995l )), who calcu- 
lated the the r mal eq uilibrium gas properties of a diffuse ISM. 
IWolfire et all i|l995h found that a stable two phase medium 
was produced, with a transition from hot material at den- 
sities < 10 -0,5 cm~ 3 to cold, molecular material at densities 
> 10 1 ' 5 cm" 3 . The Jeans mass of the warm (T = 10 4 K) 
phase is approximately 1 x 1O 8 M0, whereas the Jeans mass 
of the cold molecular gas (T = 100K) is ~ 2000M©. Typical 



cosmological galaxy simulatio ns currently have mas s resolu- 
tions no better than 10 6 M Q l|Okamoto et all l|2005l )). many 
orders of magnitude away from being able to resolve the 
relevant scales for accurate tracking of the cold molecular 
phase. Note that these simulations do however resolve Mj 
in the warm phase. Given this limitation, we are forced to 
instead introduce another particle type in our simulations, 
called a 'sticky particle', to represent the clouds. These move 
as ballistic particles through the ambient medium, yet when 
they encounter another sticky particle interact inelastically 
based on a set of collision rules that mimic the coagulation of 
clouds. Stic ky particles have been used before for similar rea - 
son s by e.g. iLin fc Pringle! Jl976l) ; Ijenkins fc Binnevl (jl994T ) 
an d lSemelin fc Combes! h20oK and seem to have been intro- 
duced originally by iLarsonP 19781) and iLevinson fc Roberts! 
|l98ll ). We show below that our imposed collision rules pro- 
duce cloud statistics that are similar to those determined in 
nearby galaxies, for which they can be measured. 

We begin by giving a brief overview of the current the- 
ory of the ISM and the formation of clouds, on which our 
recipes are based (section [2}. We then introduce the physics 
we model with the sticky particle prescription (section [3} , 
and constrain all of the free parameters in the model (section 
[4| . We then present results from simple isolated galaxy sim- 
ulations (section [SJ) and investigate the effects of changing 
the physics in the sticky particle model (section 15. 3[) . Sec- 
tion [6] contains conclusions and details of future work. The 
ISM is a complicated system and the nomenclature used 
to describe it is correspondingly complex. For this reason 
appendix [A] contains a list of symbols and their associated 
meanings. Readers not interested in numerical details can 
read the summary of the model in section [3] and then di- 
rectly proceed to the analysis of the properties of the simu- 
lated galaxies presented in section [5] 



2 STAR FORMATION IN DISK GALAXIES 

2.1 The interstellar medium in disk galaxies 

According to the models of McKcc & Ostrikerl (1 19771 ) (here - 
after M077; see also lEfstathioul (|2000D ; iMonacol l|2004l) ; 
iKrumholz fc McKed l|2005h ) the ISM of the MW consists of 
at least three separate and distinct gas phases: a hot tenuous 
medium at a temperature of ~ W e 'K; cold, dense molecular 
clouds at a temperature of ^ WOK and a warm medium 
that exists at the boundaries between clouds and the hot 
medium with T ~ 10 K. In the MW, the hot medium has 
a filling factor of 70-80% and the cold clouds account for 
a few percent of the volume (M077). Different techniques 
are used to observe these different phases, with radio ob- 
servations probing roto- vibrational transitions of molecules 
(mainly CO), the 21-cm line probing atomic hydrogen, and 
UV- and X-ray obse r vation s probing the hot phase, see e.g. 
iBinnev fc Merrifieldl { 19981 ) for an overview and further ref- 
erences. The fact that different techniques are used to ob- 
serve the different gas phases might exaggerate the degree 
to which these phases are really distinct. 

Observations of star formation in the MW show that 
most stars form in groups, either as gravitationally bound 
clusters or unbound associations, in the most massive of 
the molecular clouds (Giant Molecular Clouds, hereafter 
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GMCs), with masses ~ 10 6 Mp ), and sizes of order 30- 
50pc see e.g. iBlitz fc Thaddeusl (|l98Ch and lLada fc LadH 
l|2003! ) for recent reviews. The relatively small sizes of 
GMCs currently limit detailed observations of such clouds 
to nearby gala x ies, w it h recent surveys done in t he M W 
llSolomon et all dl987h: iHever. Carpenter fc Snel] ll200lh) 
M33 jRosolowskv fc Blitj ||2004j n and the LMC (|Fukui et a] 

teooiR 



Blitz et a] l|2006l ) present detailed observations of 



GMCs in five local galaxies. GMCs are excellent tracers of 
spiral structure to the extent that they are often used to 
define the location of arms, in a similar way as, for exam- 
ple, HII regions or massive stars are. There is a good corre- 
lation between the locations of GMCs and filaments of HI, 
although HI is often observed without accompanying molec- 
ular clouds, suggesting that the clouds form from HI. 

Observations in more distant galaxies are currently lim- 
ited to measuring the mean surface density of molecular gas 
with more detailed observatio ns awaiting new instrum ents 
such as ALMA0. According to I Young fc Scovillel l|l99ll ). the 
fraction of gas in the molecular phase depends on Hubble 
type, with early-type spirals tending to be predominantly 
molecular and late-types atomic. Optically barred spirals 
show a clear enhancement of CO emission along the bar. 
This suggests that the large-scale structure of the galaxy 
affects the formation of clouds and hence, presumably, also 
the star formation rate. A physical implementation of star 
formation should therefore aim to track the formation and 
evolution of molecular clouds. But how do GMCs form? 



2.2 The formation of molecular clouds 



Altho ugh the early models by iField. Goldsmith fc Habind 
( 1969) assumed that the different phases of the ISM were in 
pressure equilibrium, modern observations paint a picture 
of a much more complex and dynamic situation in which 
the ISM is shaped by turbulence, possibly powered by SNe 
and the large-scale dynamics of the galactic disk itself, see 
e.g. lBurkertll|2006l ) for a recent review. A galaxy-wide effect 
seems to be required to explain the observed Hubble-type 
dependence of cloud properties. 

Yet how GMCs f orm in this complex enviro nment is 
not well understood l|Blitz fc Rosolowskvl l|2004 )'). Some 
authors have suggested that GMCs for m by the coag- 
ulatio n of pr ee xistin g molecular clouds l|Kwan fc Valdej 
(|l983h ; lOortJ (|l954T )l while others have argued that 
GMCs form primarily from atomic gas throu gh insta- 
bility or large-scale shocks (|Blitz fc Shul dl98d ')'). A vi- 
able mechanism by which this could occur is the forma- 
tion of convergent flows induced by a passing spira l arm 
ijBallesteros-Paredes Vazquez-Semadeni fc Sclol {l999)). Of 
course, both modes of formation may occur: in high density 
regions, where the vast majority of hydrogen is molecular it 
seems likely that molecular clouds form from the coagula- 
tion of smaller clouds, whereas in the outskirts of galaxies 
where the gas is predominantly atomic the compression of 
gas in spiral density shock waves provides a more plausible 
formation mechanism. 
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inter arm re gions, formed from atomic gas by shocks, as 



Observationally, star formation in disks se ems to oc- 
cur on ly above a given surface-density threshold (jKennicuttl 
(1989)), with star formation dropping abruptly below the 
threshold even th o ugh t he gaseous disk may extend far 
beyond it. ISchavei (2004) describes a model in which this 
threshold arises naturally due to the thermal instability 
when gas cools from 10 4 K to the cold phase (~ 100K), 
rendering the disk gravitationally unstable on a range of 
scales. This suggests that a combination of thermal instabil- 
ity and large-scale dynamics may be responsible for deter- 
min ing when and whe re GMCs form. 

lElmegreenl ([2000) discusses observational evidence that 
GMCs are in fact short lived entities that form, make stars 
and disperse again on their dynamical time scale. This 
short time-scale alleviates the need for an internal energy 
source to sustain the observed internal supersonic turbu- 
lence, somstlriiig_tli^^ for a long 
time. iPringle. Allen fc Lubowl (|200ll ) discuss this assump- 
tion in more detail and suggest that GMCs form from ag- 
glomeration of the dense phase of the ISM, already in molec- 
ular form, when compressed in a spiral shock. They envisage 
the pre-existing molecular gas to be in dense 'wisps', in the 
eg 

simulated bv lKovama fc Inutsukal (|2000h . 

Recent numerical simulations support the view that 
when clumpy gas is overrun by a density wave it produces 
struct ures that resemble GMCs. IWada. Meurer fc Normar] 
(2002) present high-resolution two-dimensional simulations 
of the evolution of perturbations in a cooling, self-gravitating 
disk in differential rotation. They show that the disk devel- 
ops stationary turbu lenc e, even without any stella r feedback. 
iBonnell et aj i|2006l ) and lDobbs fc Bonnelll (|2006h performed 
three dimensional simulations of the passage of clumpy cold 
gas through a spiral shock. Their simulations produce dense 
clouds, with large internal velocity dispersion, reminiscent 
of the 'supersonic turbulence' seen in GMCs. They note 
that the velocity dispersion is generated on all scales si- 
multaneously, in contrast to what is usually meant by tur- 
bulence where energy casc ades from large to small scales. 
iMac Low fc Klessenl (|2004h review the current state of the 
art in simulations and models of GMC formation, including 
references to more recent work. In this picture of clouds, 
GMCs are temporary structures formed and dissolving in 
converging flows. They do not require an internal source of 
energy, are not in virial or pressure equilibrium and need 
not even be gravitationally bound. They point out that the 
relative contribution of galactic rotation and stellar sources 
to driving the observed turbulence is not clear. 

2.3 Star formation in molecular clouds 

In the older picture of cloud formation, GMCs were long- 
lived, gravitationally bound, virialised objects. The presence 
of supersonic turbulence ensures that clouds do not imme- 
diately collapse to form stars, as this would predict a star 
formation rate for the MW which is far higher than ob- 
served. Locally unstable clumps collapse to form proto-stars, 
which built-up their mass to produce the initi al mass func- 
tion (I MF) through competitive accretion (e.g. lBonnell et al 

fiaaa))- 

However, simulations show that the energy contained 
in supersonic motions is quickly dissipated even in the pres- 
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enceof magnetic fields (see references in lMac Low fc Klessenl 
l|2004 )). To sustain the turbulence therefore requires an 
energy source, for example from star formation, yet some 
clouds have turbulence but no current star formation. 

The modern picture is one in which clouds are short- 
lived structures and the turbulence results from the same 
process that formed the cloud in the first place. Observa- 
tionally, GMCs turn a small fraction e* ~ 0.1 of their mass 
into stars before they disperse again. This low star forma- 
tion efficiency of clouds may be due to the fact that they are 
short lived. The short life times of (star forming) clouds also 
follows from the sm all age spread in star clusters (see e.g. 
I Gomez et all (fl992)). and indicates that star formation in a 
given cloud only lasts for a few million years. The short life- 
time s of GMCs then also suggests that competitive accretion 
(e.g. iBonnell et all (Il997l)\ is less important in shaping the 
IMF jPadoan fc Nordlundl (120021 ) 1. 

Turbule nce generates a rang e of su bstructures inside 
a GMC, and IPadoan fc Nordlundl (|2002h suggest that such 
'turbulent fragmentation' builds a mass spectrum of proto- 
cores, some of which will collapse under their own gravity to 
form stars. The resulting IMF is a power-law due to the self- 
similar nature of the turbulence. Only cores dense enough so 
that self gravity can overcome their magnetic and thermal 
energy can collapse. This consideration flattens that IMF at 
low masses, and prevents very low-mass cores from forming 
stars. They also argue that the maximum mass is a fraction 
of the overall cloud mass. 

A young stellar population does, of course, dump a lot 
of energy into its surroundings through stellar winds, ion- 
isation, and SN explosions. Even if these do not drive the 
observed turbulence, they may contribute to the destruction 
of the cloud, and prevent further star formation. Most sim- 
ulations use such feedback from star formation to regulate 
the star formation rate. 



2.4 Summary 

The current theory of star formation in disk galaxies sug- 
gests that supersonic turbulence, generated by a combina- 
tion of galactic rotation and SNe, regulates the formation 
of proto-stellar cores inside more massive molecular clouds. 
Some fraction of these clouds can form stars, before the 
cloud itself is destroyed, by a combination of stellar feedback 
and the turbulence that built the cloud in the first place. 
The cooling of the cores to temperatures ^ 100K is domi- 
nated by grains and CII fine-structures lines and is opposed 
by photo-electric heating from sma ll grains and poly cy clic 
aromatic hydrocarbons (PAHs), see lWolfire et all (|l995l ). 

A numerical implementation of these processes requires 
high resolution to model the interstellar turbulence and fol- 
low the contraction of cores of masses a few Mq at densities 
above 1 particle per cm -3 . Such challenging simulations may 
be feasible in the near future for high-redshift galaxies but 
are not possible yet for z — galaxies. Below we describe a 
model of cloud formation that tries to incorporate some of 
these processes with some simple rules. 

We would like to apply the same rules to cosmologi- 
cal simulations of galaxy formation, which clearly requires 
a leap of faith. High-redshift galaxies may not have a well- 
defined disk, and hence the properties of the supersonic tur- 
bulence and the GMCs may well be very different. In the 



high-redshift simulations of lAbel. Bryan fc Normanl (2000), 
the first molecular cloud is close to hydrostatic equilib- 
rium, with pressure support slowly leaking away as it cools 
through molecular hydrogen line emission. This quasi-static 
evolution, reminiscent of a cooling flow, is very different from 
the dynamic turbulent fragmentation envisaged in the MW, 
with corresponding large differences in the predicted IMF. 
Furthermore, if the properties of the GMCs were similar, 
the behaviour of the cores may still be very different, with 
the reduced grain and metal cooling at higher z, making for 
a different IMF. 

Despite these difficulties we must start somewhere. The 
assumption that the physics of the ISM (and therefore the 
stellar IMF) is similar at redshift zero and in the high 
redshift universe is common in the simulation community. 
The understanding gained through these necessarily simpli- 
fied simulations will allow us, over time, to investigate the 
physics relevant to galaxy formation at higher redshift. 



3 DETAILS OF THE MODEL 

As demonstrated in section [TJ in typical simulations of 
galaxy formation we can resolve the Jeans length of the am- 
bient gas phase and so treat its hydrodynamic properties 
consistently. However, we cannot yet resolve the properties 
of the cold molecular phase of the ISM. We therefore follow 
the evolution of the ambient gas phase using a hydrodynamic 
simulation code, whereas we treat the cold phase using a sta- 
tistical model that encapsulates the physics relevant to the 
formation and evolution of molecular clouds. In this section 
we introduce the properties of the sticky particle model and 
describe the p hysics we have im plemented. 

Following lEfstathioul (|2000l ) we consider the ISM to con- 
sist of warm and hot ambient materials, and cold molecu- 
lar clouds. We additionally treat the properties of SN rem- 
nants. Throughout this paper the properties of the ambient 
medium will be represented with the subscript h, the prop- 
erties of the molecular clouds with the subscript c, and the 
properties of the gas internal to SN remnants, or hot bub- 
bles, with the subscript b. 

The ambient gas phase is represented using th e entropy 
conser v ing, parallel Tree-SPH code GAD GET2 (|Springell 
l|2005l ); ISpringel. Yoshida fc White] l|200ll )), which is a La- 
grangian code used to calculate gravitational and hydro- 
dynamic forces on a particle by particle basis. Smoothed- 
P a rticle Hydro dyna mics (SPH) was originally in troduced 
by iLucvl (Tl977l) and iGingold fc Monaghanl (119771 1. see e.g. 
iMonaghanl (|l992l 1 for a review. We will refer to the gas com- 
ponent treated using SPH as ambient gas to distinguish it 
from the cold molecular gas. Ambient gas at temperatures 
around W 4 K will be referred to as 'warm' , and gas at tem- 
peratures of lO 6 -^ and higher will be called 'hot' . We will see 
that in galaxy formation simulations, this ambient (i.e. non- 
molecular) medium naturally develops three relatively well- 
defined phases: a warm (T ~ 10 4 K) component in a galactic 
disk, a hot (T ~ 10 6 K)) tenuous component of shock-heated 
gas in the halo, and a similarly hot component resulting from 
gas heated by SN. The fourth, cold (T ~ 100K) and molec- 
ular cloud phase is represented with sticky particles, which 
interact gravitationally with all other material in the simula- 
tion and are allowed to stick together forming more massive 
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Figure 1. Summary of the physical processes that operate in our 
model of a two phase interstellar medium. The boundaries be- 
tween molecular and giant molecular clouds and between heated 
and non-heated diffuse gas are somewhat arbitrary and they are 
separate out in this figure only to highlight the different physical 
mechanisms that are operating at any given time. 



sticky particles. Stars and dark matter are both treated as 
collisionless particles by GADGET2. 

The different phases of the ISM may interact with each 
other as follows: thermally unstable ambient gas may col- 
lapse into molecular clouds via thermal instability (section 
13 - 1 1) - Molecular clouds can interact with each other to form 
GMCs ("section l3.2p . GMCs then collapse into stars (section 
I3.3f) . Stars disrupt the cloud they formed from and may, via 
SN feedback, return energy (section I3.4.1f) to the ambient 
phase. Hot bubbles blown by SNe can evaporate cold clouds 
(|3.5|) and heat the ambient medium. 

Fig [T] contains a summary of all of the physics imple- 
mented in our model. Arrows represent a transfer of mass 
and/or energy from one phase to another. The distinction 
between clouds and GMCs is somewhat arbitrary; they are 
separated in the figure to allow an easy pictorial represen- 
tation of mass transfer within a single phase. Appendix IA1 
contains a list of frequently used symbols and their meaning. 

Each physical process will be treated in turn in the re- 
mainder of this section. We first introduce the physics rele- 
vant to each physical process before discussing the numerical 
implementation. We will also give our preferred physical val- 
ues for the various parameters that occur. How we choose 
these is discussed in section [4] 



3.1 Radiative Cooling And The Formation of 
Molecular Clouds 



Begel man fc McKeei (|l990h show that under appropriate 
physical conditions, a thermal instability may operate in the 
ambient gas, which causes a fraction of the gas to condense 
into much colder molecular clouds. The sticky particle star 
formation prescription contains a basic representation of this 
process, based on a detailed treatment of baryonic radiative 
cooling. 



3.1.1 Relevant Physics 

The radiative processes that we take into account are 
Compton cooling off the microwave background, thermal 
Bremsstrahlung cooling, line cooling and photo-ionization 
heating from Hydrogen, Helium and metal species in the 
presence of an imposed ionising background. These routines 
were developed for a different project and will be described 
elsewherfl Briefly, they use tabulated rates for radiative 
cooling and photo-ionization heating for many species and 
ionization states computed assuming ionization equilibrium 
usi ng cloudy (version 05.07 of the code last described 
bvlFerland et al.1 il 19981) ) with a UV background given by 
lHaardt fe Madaul |200lT ). The rates are tabulated element 
by element and we will assume solar abundance ratios and 
specify a fixed metallicity of the gas in solar units. We do, 
however, note that the behaviour of the system may depend 
upon precisely which value of the metallicity we choose, and 
investigate this in section IS31 

Other processes such as cosmic ray heating, and cooling 
by dust and atomic lines that affect the molecular gas in 
clouds are not treated explicitly since we do not model the 
internal properties of the clouds themselves. 

We include a simple model to determine the rate 
at which the ambient gas forms molecular clouds. When 
we identify ambien t gas that is thermally unstable 
l|Begelman fc McKeei l|l990l )) we allow it to collapse into 
molecular clouds. The rate at which this process occurs is 
goverened by the rate at which the gas is losing thermal 
energy by radiative cooling. 



3.1.2 Numerical Implementation 

Following lYepes et all (Il997l ) we define a density threshold, 
p t h, to determine when gas becomes thermally unstable. Gas 
with p < p t h undergoes ordinary radiative cooling. Gas with 
density above the threshold becomes thermally unstable and 
begins to be converted to molecular clouds. In addition to 
the density criterion we add a maximum temperature (T t h) 
for gas to be called thermally unstable which has the effect 
of preventing SN heated gas in dense regions from collapsing 
straight to the cold phase. 

When gas has been identified as thermally unstable it 
begins to form molecular clouds at a rate controlled by the 
rate at which the gas can lose thermal energy by radiative 
cooling 



dt 



dp h 
dt 



1 



Uh 



-A net (p h ,u h ) 



(1) 



where u represent an internal energies per unit mass. The 
subscripts h and c refer to the ambient phase (either warm 
or hot) and cold phase respectively. A Ilct is the cooling rate 
of the ambient gas (ergs cm -3 s _1 ). We assume that the cold 
clouds remain at a fixed temperature of T c = 100K hence 
their thermal energy u c is a constant as well. 

In practice, each ambient gas particle is identified as 
either thermally unstable, or non-thermally unstable, non- 
thermally unstable gas undergoes radiative cooling; ther- 
mally unstable ambient gas forms molecular clouds at a 

2 We would like to thank our colleagues J Schaye, C Dalla Vecchia 
and R Wiersma for allowing us to use these rates. 



6 C. M. Booth, Tom Theuns and Takashi Okamoto 



rate controlled by the radiative cooling rate, as described 
byEq.IU 

In this way each ambient gas particle can keep track 
of what fraction of its mass is in the form of molecular 
clouds. Gas in the molecular phase is ignored for the pur- 
poses of the SPH calculation. When the amount of mass in 
the molecular phase in a particle reaches the resolution limit 
of the simulation a separate 'sticky particle', representing 
many sub-resolution molecular clouds is created. This pro- 
cess decouples the molecular clouds from the associated am- 
bient phase. Since we cannot resolve the individual molecu- 
lar clouds in each sticky particle we work with the mass func- 
tion of clouds. Initially we assume that the molecular clouds 
formed through instability are all in the smallest mass bin, 
that is that the clouds formed by thermal instability are very 
small, and will interact to form more massive clouds. In the 
following section we describe the behaviour and evolution of 
the sticky particles in the simulation. 

3.2 Cloud Coagulation and GMC Formation 

Molecular clouds are typically many orders of magnitude 
more dense than the medium they form in (M077), and 
their behaviour is governed by a different set of rules than 
the ambient medium. This section describes the physics of 
the simplified molecular clouds in the sticky particle model 
and how it is implemented. 

3.2.1 Relevant Physics 

We assume that clouds may be treated as approximately 
spherical objects that obey a power law relation between 
mass (M c ) and radius (r c ) 

/ Mr \ °- 3 

= 36 (wk) pc - (2) 

Here, a c describes how clouds grow as mass is added to 
them (if they remain at constant density then a c — 1/3), 
and M re f and r ro f are a reference mass and radius used to 
fix the normalisation of this relation. The lower bound on 
molecula r cloud masses is typically calculated to be WOMq 
l|Monacd (|2004r i) due to the efficient destruction of smaller 
molecular clouds by photoionization. We introduce an upper 
limit by converting molecular clouds with large masses into 
stars (see section 13.31 for discussion) . In order to facilitate 
easy estimates of the relative importance of various effects 
we have substituted typical numbers and units into most of 
the equations in this section. 

3.2.2 Numerical Implementation 

Each sticky particle represents numerous cold clouds. Sticky 
particles are hydrodynamically decoupled from the ambient 
SPH phase of the gas and interact only gravitationally with 
the other phases in the simulation. However, when two sticky 
particles collide they may coagulate to form a more massive 
sticky particle. The mass of the smallest molecular clouds 
is typically orders of magnitude below the mass resolution 
in a cosmological simulation. We represent an entire mass 



spectrum of clouds statistically inside of each sticky particle. 
Our formalism to treat the evolution of the mass function of 
clouds internal to each of the 'multiple cloud' particles will 
start from the Smo l uchow ski equation of kinetic aggrega- 
tion (|Smoluchowskil which describes the behaviour 
of a system consisting of ballistic particles that can inter- 
act via mergers. The coagulation behaviour of this system is 
driven by a coagulation kernel, K(m\,m2), that represents 
the formation rate of clouds of masses m — mi + m,2, 

K = <E« app }„, (3) 

where u app is the relative velocity of the clouds and E is 
the collision cross section. For a Maxwellian distribution 
of velocities wit h three-dimen s ional dispersion a we obtain 
(^app) = 1.3cr (|Lee fc Nelson! (l98l)). The product of the 
approach velocity and the collision cross section is averaged 
over the distribution of relative velocities. The cross section 
is 

E «,(,. + I tf(, + a0 J£±«LJ-), ,4, 

v telle Ca PP 7 

where the first term represents the collision geometric cross 
section and the second t erm r epresents the effect of gravita- 
tional focusing (Saslaw (1985)). The focused term becomes 
significant when the approach velocity is not much larger 
than the internal velocity dispersion of the system. In most 
cases of interest the geometric term will dominate so the 
focused term is neglected. In these calculations we need to 
assume that molecular clouds, although transient and turbu- 
lent, are stable for long enough for coagulation to take place. 
This is reasonable because the cloud velocity dispersion is 
typically l arger than the sound speed of the cold cloud gas 
(Monaco! l |2004 )). 

To model the cooling of sub-resolution molecular clouds 
via gravitational interaction it has been assumed that when 
molecular clouds with relative velocities, u app greater than 
"stick (a parameter in our simulations) collide they do not 
merge, but rather bounce back with relative velocity a frac- 
tion, rj, of the initial approach velocity. Clouds with relative 
velocities less than v st ick merge. For simplicity it has been as- 
sumed that the velocity distribution of clouds is Maxwellian 
with a velocity dispersion that is a function of cloud mass, 
a = a(m). 

The upper and lower bounds on the molecular cloud 
mass function are set such that the smallest mass bin is com- 
parable with the smallest clouds was can observe, and the 
largest molecular clouds are approximately the same mass as 
the largest clouds in the MW. The mass function is discrete. 
All clouds are assumed to form at the lowest mass, M m i n , 
and then the mass of each bin is a multiple of this value. This 
discrete mass function is neccessary when working with the 
Smoluchowski equation. 

In order for us to be able to hold a mass function with 
a large number of bins internal to every single sticky par- 
ticle without the requirement to store one number for each 
mass bin we parameterize the mass function as a third or- 
der polynomial, and store only the four coefficients between 
timesteps. 

As these sub-resolution clouds interact and merge, the 
one-dimensional velocity dispersion a(m) changes, which af- 
fects the rate of evolution of the cloud mass function, n(m). 
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Let E m = 3/2m<7 2 (ra) denote the random kinetic energy of 
clouds with mass m. The one-dimensional velocity disper- 
sion is related to the three dimensional velocity dispersion 
by <tid = (J'ad /v3. Em and may change due to three distinct 
processes: 

• Clouds with masses m' and m — m' merge to form extra 
clouds of mass m, increasing _B m at a rate -E ga m 

• Clouds with masses m may merge with clouds of any 
other mass decreasing the number of clouds of mass m. This 
process decreases E m at a rate -Ei oss 

• Clouds with mass m may interact collisionally with 
clouds of any other mass and so lose kinetic energy. This 
process decreases E m at a rate E coo \ 

The net change in kinetic energy for particles of mass 
m during some timestep At is given by 

AE m = At — ^-Egain — -Eloss — -Ecoolj At . (5) 

And for this change in kinetic energy, the corresponding 
change in velocity dispersion is given by 



Details of the equations used to model these processes are 
given in Appendix [B] and the method by which they are 
solved numerically in Appendix [Cl 

The same processes (cooling and merging) are followed 
explicitly for the individual particles in our simulations, 
which can interact in the same two ways as the unresolved 
sub resolution clouds by merging to form more massive 
sticky particles, or cooling to decrease their relative ve- 
locities. Following the same rules should allow us to re- 
move much of the resolution dependence of the star for- 
mation prescription. As the mass resolution of a simula- 
tion is degraded, more massive clouds will be treated with 
the sub-grid physics; our implementation should ensure that 
the large scale results are approximately the same. This is 
demonstrated in section 14.1.21 

3.3 Cloud Collapse and Star Formation 

The vast majority of stars form in Giant Molecular Clouds. 
This process is described in the sticky particle model by 
allowing the most massive clouds in the galaxy to collapse 
into stars. 

3.3.1 Relevant Physics 

We follow the process of star formation in our simulations 
by waiting for star forming clouds to be created by the co- 
agulation process described in section 13.21 We define star 
forming clouds to be clouds of a mass similar to the most 
massive clouds observed in the MW (~ 1O 6 M0). When one 
of these star forming clouds is created it is assumed to col- 
lapse on a short timescale and approximately ~ 10% of 
its mass is converted into stars, whilst the remainder is dis- 
rupted by stellar feedback processes including stellar winds, 
SN feedback and photoionization. This process reflects that 
although stars may form in less massive molecular clouds, 
it is not until the relatively rare, m assive O and B sta rs are 
created that the cloud is destroyed l|Elmegreenl l|l983l )). 



We assume that each cloud collapse for ms a sing l e stel - 
lar population with an IMF of the standard ISalpeterl (|l955f ) 
form 

N(M) dM oc M~ (1+x) dM , (7) 

where x is the slope of the IMF and takes the usual value 
of 1.35. The masses of stars are assumed to lie between well 
defined minimum and maximum values, M*, m j n and A/*, max . 



3.3.2 Numerical Implementation 

The treatment of star formation adopted in most simulations 
is to identify gas that is likely to be star-forming and impose 
a star formation rate given by the Schmidt law, 

P* = Cp^ . (8) 

Here, p* and p gas denote the rate of star formation per unit 
volume and the gas density respectively. This power law rela- 
tion between star formation rate (SFR) and gas density was 
found to hold over many orders of magnitude by iKennicuttl 
( 1998), who constrained the exponent to be Nsf = 1.4±0.2. 

We take a different approach: unstable molecular clouds 
are identified in the simulations as any cloud with a mass 
greater than M B f. We identify the formation of these massive 
clouds by using the cloud mass function, as stored internally 
to every single sticky particle. These unstable clouds are 
assumed to collapse on a very short timescale, forming stars. 

As soon as a cloud of mass M s f forms, it is assu med to be 
disrup ted by OB stars on a timescale of ~10Myr (jMatznerl 
(2002)), the rest of the massive cloud is broken down into 
smaller clouds and the coagulation process begins all over 
again as described in section 13.21 This process is modelled 
by taking the fraction of the cloud's mass that does not 
turn into stars, 1 — e», and assuming that the net effect of 
the stellar feedback processes is to fragment the CMC into 
the smallest clouds represented in the sticky particle internal 
mass function. This has the net effect of steepening the cloud 
mass function. 

Each star particle formation event represents the for- 
mation of a single stellar population of stars that are all 
assumed to have the same age, and to be drawn from the 
Salpeter IMF. Each stellar particle is therefore formed with 
a mass approximately equal to e* times the mass of a star- 
forming cloud. If this particle mass is not allowed by the 
mass resolution of a given simulation then we either store 
up unresolved stars internal to a sticky particle (if the star 
mass is too small to be allowed), or split it into multiple, 
equal mass particles (if the star mass is too large to be al- 
lowed) . 



3.4 Supernova Feedback 

Our simulations include only energy feedback from type II 
SN. These events return energy from the stars to the ambi- 
ent phase. We note that it is not currently computationally 
feasible to resolve the properties of SN remnants so we treat 
them with a simple, analytic prescription. The mechanism 
by which SN feedback is implemented in our model is dis- 
cussed here. 
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0.003 . 



(9) 



3.4-1 Relevant Physics 

Each star of mass greater than 8Mq releases 10 51 i?5i ergs in 
thermal energy when it undergoes an SN event. The lifetime, 
t, of a star of mass M (whe re M > 6.6Mq) is given by 
l|Padovani fc Matteuccil (|l993T )l 

Gyr ~~ ■ VJWq ) + 

Each SN explosion can be approximated as the injec- 
tion of energy at a single point in space. If we assume that 
the ambient density on scales of interest is approximately 
homogeneous, with density ph, the n each SN exp losion can 
be modelled as a Sedov blast wave <|Sedovl <|l959f) ). Accord- 
ing to this solution, if at time t — we release an amount 
of energy Eb, then after time t the resulting blast wave will 
have reached a radius rt given by 



(-) 



1/5 



2/5 



/^/lO^ergs 
292 U/0.1cm-3 



1/5 



(t/10Myr) 2/5 



pc. 



(10) 



These hot SN bubbles have two main effects. Firstly, as 
they expand and decelerate the SN heated gas will get mixed 
in with the surrounding ambient medium; the net result of 
this process is the heating of the ambient medium. Secondly, 
as discussed in section 1331 any cold clouds caught inside an 
SN bubble will undergo evaporation. 

There are two main assumptions that must hold for 
the Sedov solution to be valid, the pressure of the ambi- 
ent medium, and the cooling rate inside the bubble, must 
both be negligible. Often at least one these assumptions is 
invalid. If the ambient medium has a low density and is very 
hot, for example due to a previous set of explosions, then 
its pressure is no longer negligible and the Sedov solution 
breaks down. If the ambient medium is dense then radiative 
cooling becomes an important process. In the remainder of 
this section we describe various modifications to the stan- 
dard Sedov solutions, which allow us to model SN remnants 
in a wider variety of conditions. 

In the case of a hot, tenuous medium the r adius of each 
blast wave is increased (|Tang fc Wand |2005l )V These au- 
thors derive a fitting formula for the velocity of a SN blast 
in a hot medium, which is accurate to within 3% 



nit) = 



3/5 



dt' 



= 156 



i/Myr 



v + 1 



3/5 



dt' pc, 



(11) 



(12) 



where Ch is the sound speed of the ambient medium. We as- 
sumed a temperature of Tu = 10 6 K, mean molecular weight 
of p = 0.58, blast wave energy of 1 x 10 51 ergs and an am- 
bient density of 0.1 atoms per cm 3 in order to illustrate the 
order of magnitude of r%. t c is a characteristic time, 

1/3 



'2 A 5 E b 



0.012 



PhCl 



(e/i.i4) s 



£ 6 /10 51 erg 



(p h /0.1cm-3)(T^/10 6 K)5/2 



1/3 



Myr. 
(13) 
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Figure 2. Comparison between SPH simulation of a Sedov blast, 
the Sedov solution and the hot medium correction of Tang & 
Wang (2005). The points represent individual SPH particles, the 
dashed line is the Sedov solution and the dotted line is the blast 
wave radius as calculated with the hot medium correction. The 
initial condition had a density of 0.001 atoms per cm 3 and a tem- 
perature of 10 6 K. 10 51 ergs were injected to the central 32 particles 
at to- This plot was made after the blast wave had evolved for 0.3 
Myr. 



where £ equals 1.14 for a gas with adiabatic index 7 = 5/3. 



This solution matches the standard Sedov evolution, r& oc 
t 2 ^ 5 , closely until t ~ t c , after which the shell's velocity be- 
comes constant, r& oc t. This modification allows us to take 
into account that a the majority (^90%), of SNII happen 
in preheated SN bubbles (|Higdon et all (|l998t )) and, there- 
fore, the approximation that the pressure of the ambient 
medium is negligible is often incorrect. Fig ((2} shows the 
difference between an adiabatic gas SPH simulation of a SN 
induced shock- wave, the pure Sedov solution and the blast 
wave radius as predicted by the hot medium-modified Sedov 
solution from Tang & Wang (2005). 

Situations where radiative cooling are important 
may be taken into account using the prescription of 
iThornton et all (|l998h . whose high resolution simulations of 
SN explosions expanding in an ambient medium with tem- 
perature Th = 10 3 K, provide the total thermal energy in SN 
bubbles as a function of time, ambient density and metallic- 
ity. We per form bilinear int e rpola tion on the results in tables 
2 and 4 of IThornton et all (|l998h to obtain the SN bubble 
radius and thermal energy at any given time. 

Neither of these solutions treats the more general case 
of SN remnant expansion in a porous ISM, which may have 
regions of both high and low ambient density, and so we are 
not able to include the effects of SNe fully self-consistent. 
In most of our simulations we use the simple Sedov solution 
for the evolution of the SN blast waves, but note that the 
details of our prescription are uncertain. In section 15.31 we 
investigate the effects of using different implementations of 
the blast wave implementation to estimate how important 
the details of the behaviour of SN remnants are to the overall 
properties of the galaxy. Both the radiative cooling and blast 
wave velocity physics are varied. 
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3.4-2 Numerical Implementation 

By assuming that each stellar particle in the simulation rep- 
resents an entire population with the same age we can cal- 
culate the minimum and maximum masses of stars that un- 
dergo supernova events over any given time period using 
Eq Each of these supernovae is assumed to go off in a 
neighbouring gas particle (i.e. one for which the distance, r, 
to the star is smaller than its smoothing length, h, in the 
SPH formalism). We chose this particle randomly from the 
neighbours, with a weight computed from the solid angle, 
fl, it subtends on the sky as seen from the position of the 
star particle, 

fi = 4tt ( 1 ) . (14) 

V ^f(r 2 + h 2 )) 

This weighting forces that nearby hot, diffuse gas (which 
tends to have larger h, hence larger weight) is heated more 
frequently than cooler, denser parts of the ambient medium 
(which are dense, hence have smaller h). 

We do not transfer all SN energy to gas particles 
each timestep. Assuming that SN explosions are distributed 
evenly in time and space we can calculate for every ambient 
gas particle a 'porosity' of SN bubbles, Q = Vb/Va- For the 
volume associated with a gas particle we use Va = (47r/3) h 3 , 
and Vb = (47r/3) rt is the total volume of all the SN bub- 
bles in this particle. When Q is greater than a critical value, 
Qcrit ~ 1, the ambient phase is heated, else the available SN 
energy is carried over to the next time step. This ensures 
that the ambient phase is only heated when hot supernova 
bubbles make up a significant fraction of the volume. There 
are two motivations for this, firstly a given SPH particle 
cannot represent more than one phase at a given time. Sec- 
ondly simulations usually do not limit the timestep to be a 
fraction of the cooling time. Consider a warm, T ~ 10 4 K, 
SPH particle in the disk. If a small amount of SN energy is 
injected into this dense particle, it will cool very efficiently 
since the cooling rate is very high. It is only when the particle 
is heated to T > 10 6 K that the reduced cooling may affect 
the particle dynamically, so that it will move into lower den- 
sity gas, further decreasing its cooling rate, and becoming 
part of the hot, tenuous gas. Storing the available heating 
until the SN bubbles fill a significant fraction of the parti- 
cle is a way of easing the transition from warm to hot and 
makes the outcome less dependent on the timestep. 

To determine the porosity Q, we need to know the cur- 
rent radii, r&, of SN bubbles. The radius rb depends on the 
ambient gas properties and also on the available energy, Eb, 
as discussed in section [3.4. II Typically a single stellar parti- 
cle will undergo multiple SN events over a single timestep. 
Using Eqs. ((9]l and (|10fl and obtaining the SPH estimate of 
the ambient gas density at the position of the star particle 
we can estimate the average radius of all supernova bubbles 
blown by a given star particle at any time. Working under 
the assumption that the porosity of the ISM is low we cal- 
culate the radiative loss from each bubble separately. When 
the porosity of the ISM becomes Q > Qcrit ~ 1, the SN bub- 
bles are overlapping significantly and all coherent structure 
is assumed to be wiped out. The ambient gas particles are 
heated by the remaining thermal energy in the supernova 
bubbles and they are considered to disperse. The porosity is 
set back to zero. Note that using the Sedov solution implies 
we neglect radiative cooling in the remnants to determine 



the porosity, Q. However to determine how much energy is 
in the bubbl es once we decid e to he at the particle, we do use 
the tables of lThornton et al fl998) to account for radiative 
cooling in the SN shells. We believe that even though this 
treatment is not fully consistent, it does capture the main 
physics. 

3.5 Thermal Conduction 

Thermal conduction between the ambient and cold gas in the 
simulation is an important ingredient in the self-regulation 
of the star formation rate in our model of the ISM. 

3.5.1 Relevant Physics 

Thermal conduction has two primary effects. The first is 
to smooth out the temperature and density profiles inside 
SN remnants. In the strong explosion solution of Sedov, 
where thermal conduction is neglected, the temperature of 
the blast wave increases sharply towards the centre of the 
blast. This is due to the fact that the gas near the origin was 
heated by a stronger shock than that at the edges and there- 
after evolves adiabatically. The effect of thermal conduction 
is to efficiently transport heat from the centre of the blast 
to the outer cool regions. The temperature of the interior 
of the supernova blast, TJ,, is then approximatel y constant 
and e qual to the mean temperature of the blast l)Chevalierl 
l| 19751 ); M077): 

(^-) - 1 2 (^Y 3 ( n " r*r Eb ) 

\10 S KJ llOpJ l 0.1cm- 3j 4o 51 erg j ' { ' 

where n;, and Tj are the mean density and temperature in- 
side the bubble, respectively. We assume rb to be described 
by Sedov's self-similar solution. The density rib is also ap- 
proximately constant and is given in terms of the ambient 
density, rih, as 




The dimensionless number E con represents the effectiveness 
of evaporation, 

E con = ^(^) 2 / c TV- 1 , (18) 
3 Vpc/ 

jMcKeei (| 19771 ) ), and depends on a con = ft/ch (the ratio 
of the velocity of the supernova blast wave to the sound 
speed of the medium), the cloud's radius, r c , the volume 
filling factor of the cold clouds, / c i, and the efficiency of 
thermal conduction, (j> (see M077 for details). For a pure 
Sedov blast wave a con = 1.68. The presence of magnetic 
fields and turbulence may decrease <f) below its maximum 
value of <j} — 1- We compute f c i for each sticky particle from 
its current cloud mass spectrum given the assumed cloud 
mass-radius relation, Eq. ([2}. 

The second effect of ther mal co n duction is to e v apora te 
cold clouds. According to rtMcKeel i|l977l ); ICowiel l| 19771 )). 
the evaporation rate is well described by: 
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3.5.2 Numerical Implementation 

Since we store the mass function of molecular clouds inter- 
nal to each sticky particle explicitly (Sect. 13^2)1 . we can apply 
Eq. (|19jl along with Eq. @ to each cloud mass bin to cal- 
culate the total mass loss of a cloud over one timestep. The 
evaporation rate of the cloud depends on the temperature 
of the ambient gas, which is represented with SPH particles. 
However, as we discussed above, some fraction Q of the vol- 
ume of each SPH particle may be filled by hot SN bubbles, 
in which the evaporation rate of clouds may be much higher. 
Since we have computed Q, we can take this important effect 
into account. 

Consider a single molecular cloud in thermal contact 
with an ambient medium of (constant) temperature T. The 
mass of a cloud at the end of a timestep (Mf) is related to 
its mass at the start of the timestep (Mi) by: 



M f = \m}~ 



(1 



0.44T 



5/2 



r -^AtV /(1 - 



(20) 



where T is in units of 10 6 K, masses are in Mq, lengths are 
in pc and times are in Myr. 

Eq (|20p represents the mass loss rate for a single cloud 
in contact with a medium of temperature T. More generally 
in a porous medium a single cloud of mass m has a mean 
mass loss rate described by: 



M cloud = -QM bubble - (1 - Q)M, 



ambient ■> 



(21) 



where rh bubb i and m ambicnt represent the rate of mass loss 
for a cloud inside a supernova bubble and situated in the 
ambient medium respectively. 

Eq. (|2Up can be applied directly to the evaporation of 
a cloud in the local ambient medium (m amb i ont ). However 
to apply the same formula to the evaporation of clouds in- 
side of supernova bubbles we need to account for the fact 
that although the temperature inside the bubbles remains 
uniform, due to conduction, it is not constant in time, but 
decreases as the bubble expands. We therefore make the 
additional assumption that the mean temperature of the 
supernova remnant is constant over a timestep (a good ap- 
proximation after a short (~ 20yr) transient phase). Under 
this assumption Eq. (|20[) can be applied successfully to the 
more general case of evaporation in a porous medium. Eq 
Q19p and Eq © are used to show that the total mass loss 
rate for clouds of mass m in a volume Va is given by 
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(22) 



Under the assumption that Tb, the mean temperature of su- 
pernova remnants, and T a , the mean ambient temperature, 
are constant over any single timestep we can write 



/ M c \ = x (Mc_y 



(23) 



In order to calculate the constant of proportionality, A, we 
use an estimate of the mean temperature and density inside 
of a supernova remnant. These estimates were obtained by 
noting that by definition Q = Vb/Va- (Vb and Va represent 
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Figure 3. Evolution of a population of molecular clouds as they 
are evaporated by a hot ambient medium. The initial cloud mass 
function is a power law. The temperature of the ambient medium 
is assumed to be 10 K, the porosity of the medium is assumed 
to remain constant at 0.2, and the temperature of the supernova 
remnants is 10 6 K. Thermal conduction acts to preferentially 
destroy the smaller clouds. 



the total volume in bubbles and the ambient phase respec- 
tively). The mean radius of a supernova remnant is then 



Tb 



I 3QV A \i/a 
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(24) 



where Nsn is the total number of supernova explosions that 
have affected the local ambient medium (Calculated from 
equations [9] and [7| . The mean density inside the supernova 
remnants, rib, may then be calculated from Eq (|16[l and 
Eq (|17[) and the mean temperature from Eq (|15|) . 

Over a period of time At a cloud with mass Mi will 
evaporate to a mass of Mf, given by: 



M f 



(Mf 



\At 



l/(l-a c ) 



(25) 



Thermal conduction efficiently destroys smaller clouds, 
but its effects are far less dramatic on larger clouds. Fig <(3j 
shows the evolution of an initially power law mass spectrum 
of clouds in a hot medium. The energy used to evaporate a 
mass Mf — Mi of cold clouds is removed from the supernova 
remnants. 



3.6 Mass Resolution Limits 

The sticky particle model allows particles of all types to 
change their mass via processes including merging, thermal 
conduction and star formation. For this reason it is neces- 
sary for us to introduce numerical minimum and maximum 
masses on all particle types. We define at the initial time 
a characteristic mass resolution for our simulation, A/ C har, 
typically this is set equal to the mass of the ambient gas par- 
ticles in the initial conditions. Where more than one mass of 
ambient particles is present (for example in the model galax- 
ies discussed in section[5| we use the mass of the gas particles 
that will be forming most stars. We then define minimum 
and maximum particle masses relative to this characteristic 
mass scale. 

Ambient gas particles may have their mass decreased 
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by the formation of molecular clouds. If the total mass of a 
gas particle becomes less than 0.1M c h ar then it is converted 
entirely into a cloud particle. The ambient gas particles may 
also have their mass increased by the process of thermal 
conduction. If a gas particle becomes more massive than 
4M c har then it is not allowed to grow any more, and the 
evaporated cloud mass is given to a different particle. In 
practice this limit is rarely, if ever, reached as evaporating 
cold clouds effectively cools the ambient gas particles so they 
become inefficient at thermal conduction. 

Sticky particles may decrease their mass by star forma- 
tion and evaporation. If the mass of a sticky particle drops 
below 0.1M c har then it is either completely evaporated or 
completely converted into stars. Coagulation may drive the 
mass of a sticky particle to be very large. In practice this is 
not a real concern since when a sticky particle becomes very 
massive the rate at which its internal clouds coagulate also 
increases, causing it to form stars very rapidly. 

Stars have a maximum and minimum mass of 4M c i ar 
and 0.1Af c har. If a star forms with a mass greater than the 
maximum allowed mass it is split into a number of smaller 
star particles. A sticky particle may not form a star with a 
mass lower than the minimum allowed mass. In this even- 
tuality then the mass of the 'unresolved' stars is tracked 
internally by the sticky particle and added into the next 
star formation event until the total mass of stars formed 
reaches the resolution limit of the simulation. 

These particle mass limits keep all particle masses in the 
range 0.1A/ C h ar to 4M c h ar , which both minimises two body 
effects between very massive and very small particles and 
also prevents the formation of very many low mass particles, 
which are computationally very expensive to evolve. 



4 PARAMETER ESTIMATION 

The various physical processes in the star formation and 
feedback models each have associated with them physical 
parameters. Before we discuss the properties of our model in 
detail we discuss how its free parameters can be constrained. 

The free parameters that control the thermal instabil- 
ity and formation of the molecular clouds are p t h and T t h, 
the physical density and temperature at which thermal in- 
stability is allowe d to set in a n d rad iative cooling creates 
molecular clouds. IWolfire et all l|l995l) found that a diffuse 
ISM naturally settles into two stable phases, with a sharp 
cutoff between the ambient and molecular phases at a den- 
sity of approximately 1 atom per cm 3 . We use this as the 
value of p t h- A threshold temperature Tth = 10 5 K allows 
the gas in galaxies which cools radiatively to ~ W 4 K to 
collapse into clouds but prevents supernova heated material 
(typically at temperatures of 10 6 K) from forming molecu- 
lar clouds until it has radiated away most of its supernova 
energy. 

The properties of the molecular clouds themselves are 
contained in four parameters: r re f, M re f, and a c as defined in 
Eq @ and u c , the internal energy per unit mass of molecu- 
lar clouds. The first three values are set by comparison with 
observations of mo l ecular clouds in the nearby galaxy M33 
jWilson fc Scovilld (|l99Ch ): 



(26) 



Thus r ro f and M re f are assumed to be 36pc and 10 5 Mq re- 
spectively. This calibration (and an assumed a c of 0.3) sug- 
gest a radius of 122 ± 6pc for the largest cl o uds o bserved in 
the MW ( 6 x 10 6 M© (|Williams fc McKeel (|l997h V 

The properties of the stars and associated feedback are 
contained within four parameters: x, the slope of the IMF; 
-E51, the energy of each supernova blast in units of 10 51 erg; 
M* im i n , the minimum star mass; and M*, max , the mass of 
the largest allowed stars. For E51 we use the fiducial value 
of 1.0 noting, however, that the value of E51 is very uncer- 
tain and may be significantly higher. The effects of varying 
-E51 are investigated in section 1531 For the purposes of this 
work uncertainties in the IMF are neglected and x is as- 
sumed to_t_ake_oii_th^_^andard Salpeter value of 1.35. We 
follow iKawata fc Gibson] (120031 ) in adopting values 0.2Mq 
and 60Mq for the minimum and maximum stellar masses, 
respectively. 

The star formation efficiency in a single cloud collapse is 
also some what uncertain and i s know n to be approximately 
e* » 11% l|Williams fc McKeel l| 19971 )) in the MW. 

The thermal conduction efficiency is characterised by 
two numbers: a coll , the ratio of the blast wave velocity to 
the ambient sound speed and 4>, the efficiency of thermal 
conduction. Following M077, the value of a con is set to 2.5 
(for the ideal Sedov blast wave case, a con is 1.68, the pres- 
ence of thermal conduction changes this value). The thermal 
conduction efficiency parameter is assumed to be (f> — 1. The 
presence of magnetic fields and turbulence may change (f> sig- 
nificantly; we investigate the effect of moving away from this 
value in Sect. 15.31 

This leaves « s tick (the maximum relative cloud velocity 
for mergers) and 77 (the fraction of a cloud's velocity lost 
per non-merger collision) as free parameters that are hard 
to constrain via observation. It is noted that the large scale 
behaviour of a given simulation is largely independent of the 
value of rj. This is because the cold cloud velocity dispersion 
is always limited by w s tick- In the following section simple 
simulations are used in order to calibrate the properties of 
the physical model. 



4.1 One Zone Simulations 

4-1.1 Simulation Details 

A 'one zone model' is a periodic box that represents a fixed 
mass and volume (i.e. a static, periodic box with no mass 
outflow). The ambient ISM phase is assumed to be homo- 
geneous. Initially, for a chosen mean density of matter we 
assume that 50% of the material is initially in the hot phase 
at a temperature of To = 10 6 K. The remaining gas is ini- 
tially in cold clouds with an initial mass function that is a 
very steep (N(M)dM oc M~ s dM) power law. Numerically 
we represent the different phases as follows: the ambient 
phase is assumed to be homogeneous and isotropic and so is 
represented by a single density and temperature throughout 
the whole periodic volume, molecular clouds are represented 
by discrete sticky particles that are spawned at a random 
point in the computational volume with a random velocity, 
stars are not tracked individually, and are assumed to heat 
the whole volume evenly when they undergo SN explosions. 
The mass resolution of the molecular phase is approximately 
10 7 Mq, although this figure is varied in section 14.1.21 
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Figure 4. Star formation rate as a function of time for a one zone 
box with three different values of Ujtick- The initial gas density is 
no = 2 cm -3 in each case. Each curve follows the same general 
shape, there is an initial delay during which the first GMCs are 
forming. The unopposed collapse of the first GMCs causes a burst 
of star formation, which is quickly regulated by the effects of 
feedback from stellar winds and supernova explosions. After this 
initial burst the star formation rate in the simulation settles down 
and gradually decreases as the gas in the box is used up. 



This initial situation represents hot, dense gas that has 
just begun to experience a thermal instability and started 
forming its first molecular clouds. The volume we simu- 
late is one cubic kpc. The hot phase will evaporate cold 
clouds through thermal conduction, and can cool via ra- 
diati ve processes using a s i mple tabulated cooling function 
from ISutherland fc Dopital (ll99of ) (assuming solar metallic- 
ity). Cold cloud particles are scattered randomly through- 
out the volume and given random velocities. Clouds do not 
feel gravitational forces. Depending on the parameters, the 
ambient phase will cool radiatively to form more molecu- 
lar gas. Clouds will coagulate to form GMCs, which in turn 
form stars. The associated SNe evaporate smaller clouds, 
and may heat the ambient medium and quench the star 
formation. This sequence of events is plotted in Fig Q, 
the same general shape is observed for each value of fstick, 
there is a brief delay as the first clouds coagulate to form 
GMCs, these clouds them collapse and form stars, which 
undergo supernovae and quench the star formation. On a 
longer timescale, the quiescent star formation rate slowly de- 
creases as the available gas is consumed by stars. Since the 
dynamical equilibrium is reached on a very short timescale, 
typically a couple of hundred Myr, we assume instantaneous 
recycling when considering supernova feedback. The role of 
T t h is suppressed in the one zone simulations, due to the 
face that energy injected via supernovae cannot escape the 
volume. 

The lack of self gravity does not affect the global prop- 
erties of the volume significantly. From Eq Q assuming 
typical cloud properties (M c « 10 Mq, r « 50pc) and a 
reasonable velocity dispersion (am = 7km/s) the ratio of 
the geometric part of the cross section to the gravitationally 
focused part is approximately 0.1, therefore direct collisions 
between clouds account for the majority of the collisions and 
gravitational focusing makes for an effect of only 10%. In the 
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Figure 5. Mass function of clouds after lGyr in a one zone model. 
The dashed and dotted lines represent the slopes of the mass 
functions in the MW (Solomon et al (1987))and M33 (Rosolowsky 
& Blitz (2004)). The numbers in the legend represent the power 
law slopes in each of the galaxies. It is clear that we obtain a 
good agreement between our model and the cloud mass spectrum 
in real galaxies. 

following section we will choose a value of v st i c k by compar- 
ing the star formation rates in one zone volumes with the 
Schmidt law, and also look at the properties of one zone 
volumes. 

As noted in section|4]the properties of the simulation are 
largely independent of r\. We assume a value of 0.5 through- 
out the rest of this paper. 

4-1.2 Calibrating the base model 

The one zone model provides a useful sandbox in which we 
can investigate a wide variety of parameter choices in a rela- 
tively computationally inexpensive environment. In the fol- 
lowing section we discuss our choices for the values of the 
different parameters. The effects of moving away from this 
'base model' are discussed more fully in later sections. 

The parameters that are available for tuning the output 
of the model are as follows: the cold cloud reference size 
and radius (r re f, M re f); the slope of the cloud mass-radius 
relation, a c ; the efficiency of star formation in any given 
cloud collapse, e*; the maximum relative cloud velocity for 
merger (w s tick) an d the amount of energy ejected per SNII 
event (-E51). 

Even though the initially assumed cloud spectrum, 
N(M)dM oc M~ s dM, is very steep and far from equilib- 
rium, SNe feedback and cloud coagulation rapidly build a 
mass spectrum N(M)dM oc M~ a dM with a k 2 (Fig 
([5|), close to what is the observed cloud spectrum in the 
MW (dashed line) and M33 (dotted lines). This gives us 
confidence that, although the modelling of cloud formation 
is simple, it does produce a realistic cloud spectrum. 

The ISM model can also reproduce the observed 
Schmidt law. We find that in our model the interaction be- 
tween the coagulation of clouds and their destruction by 
stars leads to a SFR-density relation that is in good agree- 
ment with observation (Fig (|6])) if w s tick is set to 7km/s. 
This represents a reasonable value for the molecular cloud 
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Figure 6. Schmidt law. The diagonal dashed line represents the 
observed star formation law (Kennicutt 1998) and the vertical 
line represents the observed cutoff in star formation (IOMqpc - 2 ; 
Schaye 2004). Each point represents the star formation rate av- 
eraged over a period of 500Myr for a separate one zone simula- 
tion. Data is shown for three different values of tigtickj the base 
value used in all subsequent simulations is 7km/s. We calculate 
star formation rates by averaging the star formation rate in the 
simulation volume over a 500Myr period. Surface densities were 
calculated from volume densities by assuming a disk of thickness 
of lkpc. 
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Figure 8. Temperature and density of the ambient phase of a one 
zone model for a variety of different choices of initial temperature 
and density. The total gas density, ambient gas plus clouds, is 
always 2 cm -3 . The interplay of supernova feedback and radiative 
cooling quickly brings the system into an equilibrium independent 
of the initial value. 
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Figure 7. Star formation rate as a function of time for one zone 
models with three different mass resolutions. The star formation 
rate remains almost unchanged over two orders of magnitude in 
mass resolution. The coarsest mass resolution of 1O 9 M0 corre- 
sponds to the entire one-zone system being represented with a 
single particle with all clouds interactions modelled with the co- 
agulation equations. 



velocity dispersion, as considering theoretical models for the 
origin of random motion on molecular clouds, we would 
expect typical velocit y dispersions in the range 5-7km/s 
i|Jog fc Ostrikerl i| 19881 )). 

The effect of changing the mass resolution of the simu- 
lation over two orders of magnitude is demonstrated in Fig 
(J7J). Sub-resolution clouds that are simulated only by inte- 
grating the coagulation equations are designed to behave 
in exactly the same way as the resolved cloud particles in 
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Figure 9. The large scale properties of a one zone model with 
initial density no = 2 cm -3 . The physical parameters used in this 
model are the same as the base model as discussed in section [TL2] 
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Figure 10. Same as Fig 10 except with an initial gas density 
of no = 16 cm -3 

the simulation, and so we expect the simulations not to de- 
pend strongly on particle number. This is borne out by the 
good agreement between simulations carried out with only 
one resolved particle (Fig ([7]), line with mass resolution of 
10 9 Mq) where all of the physics is followed by integrating 
the sub-grid equations in a single particle and simulations 
with a hundred particles that are followed explicitly. 

As stated in previous sections, the behaviour of a one 
zone model is virtually independent of its initial tempera- 
ture and the fraction of the gas that starts off in the cold 
phase. This behaviour is demonstrated in Fig A one 

zone volume with total initial density of no = 2cm -3 was 
evolved with a variety of different initial values for the ini- 
tial temperature and initial fraction of the mass in the hot 
phase. We observe that regardless of the initial choices for 
these two quantities the system quickly settles down to its 
equilibrium state. This process occurs through the opposing 
actions of thermal conduction and supernova feedback. 

Fig (|9|) and Fig (|10[) show the behaviour of the large 
scale properties of two different one zone volumes as a func- 
tion of time. The only difference in the initial conditions of 
the two one zone volumes is their initial density Fig ([9} 
shows the evolution of a one zone volume with a total den- 
sity of 2 atoms/cm 3 ; Fig (|10p shows exactly the same plots 
for a density of 16 atoms/cm 3 . The initial temperature of 
the hot phase in both simulations is 10 6 K. In both cases the 
star formation rate follows the same general shape. There 
is a small period of time at the beginning of the simula- 
tion where small clouds are coagulating and there is no star 



formation. When GMCs are formed there is a large burst 
of star formation that is quickly quenched by feedback SN 
and thermal conduction in SN bubbles. The temperature of 
the diffuse phase is regulated by a combination of supernova 
feedback (acting to increase the temperature) and radiative 
cooling. Due to the fact that we do not allow mass to leave 
the one zone volume and also assume instantaneous recy- 
cling, the temperature profile very closely matches that of 
the star formation rate. It is noted that in the one zone 
simulation with the largest density, the temperature of the 
ambient phase is held at a higher temperature by the ac- 
tion of supernovae. The fraction of the gas in the molecular 
phase is lower in the high density simulation due to the in- 
creased amount of evaporation by thermal conduction in the 
high temperature ambient phase. In the following sections 
the star formation and feedback prescriptions are tested in 
a more realistic situation. 



5 RESULTS 

As a more physically interesting test of the star formation 
model we have conducted simulations of various isolated 
model galaxies. In this section we introduce the details of 
the two different types of simulation performed and present 
the properties of the stellar disk and associated ISM in each 
case. 

5.1 Quiescent Disk 

One of the fundamental properties that a star formation pre- 
scription must be able to reproduce is that in MW like condi- 
tions, the resulting behaviour should be similar to that in the 
MW. In this section we discuss the properties of galaxy sim- 
ulations set up to approximate the conditions in the MWs 
quiescent disk. 

5.1.1 Simulation Details 

We set up a simplified model of a MW type galaxy us- 
ing in itial conditions from GalactlCS (|Kuiiken fc Dubinskil 
(1995)). GalactlCS generates near equilibrium distributions 
of collisionless particles consisting of a disk, bulge and halo. 
These models consist of a spherical bulge component; an 
approximately exponential disk, which is rotationally sup- 
ported in the x-y plane and supported by random motion in 
the z direction; and an approximately spherical halo. 

We add baryonic material to this distribution by con- 
verting the disk and bulge in their entirety into SPH par- 
ticles at a temperature of W 4 K. 1% of the material in the 
halo is converted to baryons with a temperature of 10 6 K. 
The addition of baryonic material puts the system well out 
of equilibrium so each simulation is run adiabatically for 
50Myr to allow the galaxy to relax closer to its equilibrium 
state before the additional physics is allowed to operate. 
The total mass in the disk, bulge and halo are 1 x 1O 1O M0, 
0.43 x 10 10 M© and 1.1 x IO^Mq respectively. The mass 
resolution of particles in each of three realisations of this 
galaxy are summarised in table [1] These masses were chosen 
such that the gaseous particles in each of the three compo- 
nents have approximately the same mass and the dark mat- 
ter halo particles have masses as close as possible to the gas 
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Resolution (Mq) Disk 



Bulge 



Halo 



Low 
Base 
High 



6.6 x 10 6 5.84 x 10 6 5.6 x 10 6 
8.3 x 10 5 7.5 x 10 5 7.1 x 10 5 
1.0 x 10 5 9.2 x 10 4 1.0 x 10 s 



Table 1. Initial particle masses in three different realisations of 
the model GalactlCs galaxy that are used throughout this paper. 
The disk and bulge consist entirely of baryons, the halo addition- 
ally contains dark matter. All masses are in units of Mq. Baryons 
are added to the dark matter halo by converting a random 1% of 
the halo DM particles into gas. The dark matter particle mass in 
the halo is therefore the same as the gas particle mass. 
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Figure 11. Demonstration of the star formation properties of the 
isolated galactic disk. The continuous field represents the molecu- 
lar gas surface density of the simulated galactic disk, spiral struc- 
ture is evident. The points represent the sites of star formation 
within the previous lOMyr. Most star formation events represent 
the collapse of a single GMC, resulting in the formation of 10 5 Mq 
of stars. It is clear that star formation is occurring primarily in 
the spiral arms of the galaxy. 



particle mass. The gravitational softenings for the disk par- 
ticles (that is: gas, sticky and star) is set to O.lkpc, 0.05kpc 
and 0.02kpc in simulations GAL_LORES, GALJ3ASE and 
GALJHIRES respectively. The dark matter particles have 
softening lengths ten times larger than the disk particles. 

The GalactlCS simulations provide a test of the code in 
a situation somewhat similar to a quiescent MW disk. As dis- 
cussed in section [3] all simulations were p erformed with t he 
entropy conserving SPH code GADGET2 (jSpringell l|2005h ). 
with all of the physics discussed in section [3] implemented 
as additional modules. Table [2] contains a brief summary of 
the different simulations. Typical timestep size in one of the 
base simulations is ~ 10 4 yr, although this figure is smaller at 
early times when bursts of supernovae heat gas very strongly. 




Figure 13. A thin slice of the gas temperature field through the 
centre of a disk galaxy simulation. The arrows represent the gas 
velocity field, taking into account only gas that has been heated 
by supernovae. The generation of bipolar outflows from the galac- 
tic disk is very clear. The lower plot represents the galaxy after 
50Myr, the upper panel is the same galactic disk after 500Myr. 



5.1.2 The Base Simulations 

In this section we will discuss simulations run with the base 
set of physical parameters ( section l4.1.2p . Most simulations 
were run at the base mass resolution as defined in table [1] 
The large scale behaviour of the model galaxy is as fol- 
lows: Immediately after switching on the additional star for- 
mation physics the dense, thermally unstable gas in the disk 
and bulge collapses into cold clouds, which quickly cause a 
large burst of star formation. After approximately 500Myr 
the galaxy settles down into a quiescent state with a star 
formation rate of approximately IMq/jt. The star forma- 
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Name 


Details 


Wgas 


GALJ3ASE 


Base GalactlCs model 


19330 


GAL.LORES 


Base model with degraded mass resolution 


2450 


GAL.HIRES 


Base model with improved mass resolution 


255000 


GAL_BASE_LOSN 


E51 decreased by factor of 5 


19330 


GAL_BASE_HISN 


E51 increased by factor of 5 


19330 


GAL.BASEXOCON 


Conduction efficiency decreased by factor of 5 


19330 


GAL_BASE_HICON 


Conduction efficiency increased by factor of 5 


19330 


GAL.BASEXOZ 


Gas metallicity set to 0.5 Solar 


19330 


GAL_BASE_HIZ 


Gas metallicity set to 1.5 Solar 


19330 


ROT.BASE 


Spherical rotating collapse 


15000 


ROT.LORES 


Spherical rotating collapse 


4000 


ROT.HIRES 


Spherical rotating collapse 


45000 



Table 2. Brief table of simulation references and details. -/V gas shows the number of gas particles in the disk, bulge and halo combined. 



tion rate gradually decreases as the cold gas in the galaxy 
is consumed by stars. 

It is known that in the MW, most areas of active star 
formation are concentr ated in the galactic spiral arms (e.g. 
lEngargioIa et all (|2003l )'l. Fig (|lip shows that our simula- 
tions reproduce this behaviour. In the sticky particle model 
this behaviour occurs naturally as the converging gas flows 
in galactic spiral arms lead to an increased merger rate and, 
therefore, to the presence of more star forming clouds. Face 
and edge on temperature and density plots of the standard 
resolution galaxy are shown in Fig (|12p . The gas heated 
by SNII is preferentially situated perpendicular to the plane 
of the disk, suggesting that the feedback scheme is prefer- 
entially heating the low density gas and setting up strong 
outflows. Fig (I13|l shows the behaviour of the supernova 
heated gas in a thin slice through the centre of the disk. Ini- 
tially there is a strong burst of star formation (lower panel 
of Fig (|13l) ~). followed shortly by a burst of supernova ex- 
plosions that heat the gas around the galactic disk as hot at 
10 S K. Most of this gas is driven straight out of the halo in a 
direction perpendicular to the galactic disk. Later on, as the 
supernova rate dies down, gas is heated more gently and is 
ejected from the galactic disk in the form of a fountain rem- 
iniscent of the galactic fountains present in the MW. This 
behaviour is demonstrated in figure 1141 which shows for a 
random subset of particles from the gas disk the number of 
times they have been heated as a function of time with their 
height above the galactic disk. It is clear that upon being 
strongly heated, some particles are ejected from the galactic 
disk and fall back down a few hundred Myr later. Others 
remain in dense regions and cool immediately. Some escape 
the disk completely in the form of a galactic wind. This be- 
haviour w as also observed in the m ultiphase star formation 
models of Scannap ieco et al (|2006t ). suggesting that it is a 
more general feature of multiphase models. 

Additionally our simulated galaxies are in good agree- 
ment with the observed Schmidt law, Eq (|15fl . This be- 
haviour arises due to the self-regulation of the simulated 
ISM. At higher densities more molecular clouds are formed 
and so star formation rates are higher. Runaway star for- 
mation is prevented by various forms of stellar feedback, 
which prevent too many clouds from forming. The slope of 
the Schmidt law in the simulated galaxies may be changed 



by altering the value of M s tick' Higher values of v s tick lead to 
clouds coagulating more rapidly, and so low density regions 
of the galactic disk undergo more star formation. The effect 
is less severe in the higher density regions of the galaxy as 
strong feedback from large bursts of star formation can ef- 
fectively regulate the amount of star formation. Conversely 
a lower value of w s tick leads to a steeper Schmidt law with a 
lower overall star formation rate. This behaviour is demon- 
strated in Fig for the one zone model. The value of t? s tick 
used to reproduce the Schmidt law slope and normalisation 
in the one zone model also works in the simulated galaxy 
and the observed galaxy follows the observed Schmidt law 
very closely throughout its whole lifetime (figure (JT3J). 

The resolution independence of the star formation pre- 
scription is once again demonstrated by Fig (116[) . The star 
formation rate between the highest and lowest mass resolu- 
tion simulations is in remarkably good agreement. The gen- 
eral form followed by all simulations is that there is a strong 
burst of star formation at the initial time, this is rapidly 
quenched by supernova feedback, and a self-regulating ISM 
is set up. The star formation rate slowly increases as the gas 
in the galactic disk is either used up or ejected in the form 
of winds. 

Recent observations of the gas content of the MW have 
allowed the constru ction of maps of its gas surface density 
II Le vine et all (120061)1. In order to compare the properties of 
our model to observations, another GalactlCs model was 
generated with properties as close as possible to those of 
the MW. The total mass of the galactic disk was set to 
5 x 1O 1O M0, and the scale radius of the exponential disk to 
4.5kpc. This simulation was evolved for lGyr. The resulting 
gas distribution is shown in Fig (I17|l . O ur simulation s are i n 
good agreement with the observations of lLevine et "al (|2006r i. 

These properties suggest that the star formation and 
feedback prescriptions behave well in a quiescent disk, a 
more robust test of how they perform in a more general 
situation is given by the rotating collapse simulations. 

5.2 Rotating Collapse 

5.2.1 Simulation Details 

The second simulatio n we investigate is the co llapse of a ro- 
tating spherical halo (jNavarro fc White! (|l993h > with an ini- 
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Figure 12. A thin slice of the gas temperature and density distributions after lGyr in run GAL_BASE. The slice is taken directly 
through the centre of mass of the stellar disk. The temperature plot clearly shows regions of strongly heated gas, these are areas near to 
sites of active star formation, where the massive, short lived stars are heating the ambient material via SN explosions. 



tial 1/r density profile consisting of 90% collisionless dark 
matter and 10% baryonic material. The mass of the rotating 
sphere is 1 x 1O 12 M0 and its initial radius is lOOkpc. Veloc- 
ities are chosen such that the sphere is initially rotating as 
a solid body with a spin parameter of 0.1. Once again this 
simulation is created at three different mass resolutions, cor- 
responding to dark matter particle masses of 5.2 x 1O 8 M0, 
5.8 x 10 7 A/ Q and 2.0 x 10 7 M Q , with corresponding gravita- 
tional softenings of 1.93kpc, 0.96kpc and 0.68kpc. The rotat- 
ing collapse simulations model, in a crude way, the collapse 
of a protogalaxy, and allow us to investigate how the ISM 
model behaves when it is initially far from equilibrium and 
in situations with strong shocks and rapid density changes. 



5.2.2 The Base Simulations 

After 2Gyr the density profiles of each of the three phases 
of matter are shown in Fig (|18p and Fig (|19[) . The density 
profiles are averaged around the disk; each radial bin repre- 
sents a ring centered on the centre of mass of the disk. Figure 
(|18|) shows the radial density profiles of both the ambient 
and molecular gas. The three different resolution simulations 
once again behave in very similar ways. Fig (|19|l shows the 
radial density profile of stellar mass, demonstrating that our 
star formation prescription gives rise to an exponential disk 
and a bulge component well fitted by the standard r 1 ^ 4 law. 

The rotating collapse simulations are especially inter- 
esting because they start far from equilibrium and so fea- 
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Figure 21. Density- Temperature relation for a rotating collapse simulation showing the creation of three distinct components. In the 
final time plot region a contains gas in the halo of the galaxy. Region b is gas in the disk of the galaxy, which reaches an equilibrium 
temperature of ~ 10 4 K and c represents the approximate position of the molecular clouds, at an assumed temperature of 100K. The 
gas in the halo consists both of gas that failed to cool on to the disk, and gas that was originally in the disk but has been heated by 
supernovae. Approximately 45% of the mass is in the two hot phases, 40% in the disk and 15% in the cold molecular clouds. 



tures that arise in the final particle distribution are purely 
an effect of the physics in the simulation and not just set up 
by hand in the initial conditions. Fig (|20[l shows how the 
molecular fraction (the ratio of the mass in molecular hy- 
drogen to the mass in atomic hydrogen) varies as a function 
of distance from t he centre. Observational data from M31 
l|Dame et a l (1993)) is included as a comparison. 

The evolution of the thermal properties of the halo are 
shown in Fig (|21[) . in its initial state all the gas in the 
halo is cold. As the halo collapses it becomes dense and is 
shock heated. The gas that ends up in the disk comes to 
an equilibrium between radiative cooling and the heating 



due to SNe at approximately 10 4 K and a halo of hot, SN 
heated gas at ~ lO 6 -^ is gradually formed. The addition 
of the molecular gas (T ~ 100K, p ~ 100 - 1000cm" 3 ) 
forms an ISM with four phases: shocked halo gas, SNe heated 
material, cold molecular clouds and warm disk gas. The first 
two components are hard to di stinguish on the p — T plot. 

iBlitz fc Rosolows kv (2006) have argued that the ratio 
of molecular to atomic gas in galaxies, iZ mo i, is determined 
by hydrostatic pressure and through observations of nearby 
galaxies found the following relation: 

_ r iWfc B io.mo.QT 
Kmol ~ L (3.5 ±0.6) x 10* \ ' ( 
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Figure 14. Relationship between the number of supernova heat- 
ings and the distance from the mid-plane of the disk for a sam- 
pling of five particles from the simulated galaxy. Each different 
lincstyle represents a different SPH particle. The top panel shows 
the perpendicular distance from the centre of the stellar disk, the 
central panel shows the entropy of each particle and the lower 
panel the cumulative amount of thermal energy that has been 
dumped into the particle. It is clear that some particles with a 
higher entropy are lifted away from the galactic disk where they 
cool and rain back down on the galactic disk within a hundred 
Myr of being supernova heated. Other particles are ejected vio- 
lently from the galaxy, their density becomes very low and they 
evolve adiabatically. 



where P ex t is an estimate of the mid-plane pressure. R mo i is 
the ratio of mass in the atomic and molecular phases. Fig 
p2]l shows the data for the ROTJ3ASE simulation after 1 
Gyr alongside the observed best fit line (Eq (|27[) ). We use 
the SPH estimate of the pressure at z — 0. To calculate -R mo i 
we use all matter within a vertical distance of lkpc from the 
centre of mass of the disk and bin radially. The model of the 
ISM clearly reproduces the observed behaviour. 



5.3 Away From the Base Model 

The determination of some of the physical parameters used 
in our model is somewhat uncertain. In this section we in- 
vestigate the effect of varying some of the physics included 
in the models. Large uncertainties are present in the de- 
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Figure 15. The density-star formation rate relationship for our 
simulated quiescent disk. The diamond shaped points represent 
the observed star formation rates after 200Myr, and the triangles 
represent the star formation rates in the same disk after lGyr. 
The dashed line is the observed Schmidt law due to Kennicutt 
(1998). The points on this plot were calculated by taking radial 
bins, then calculating the mean density and star formation rate 
in each bin over a short period. The quiescent galaxy follows the 
observed Schmidt law throughout its lifetime. 
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Figure 16. Star formation rate as a function of time for isolated 
galaxy models with different mass resolutions. The mass resolu- 
tion is 64 times better in the high resolution disk than the low 
resolution one. The three simulations plotted are GAL_LORES, 
GALJ3ASE and GAL.HIRES. The fact that the star formation 
rate remains almost unchanged shows that numerical convergence 
has been achieved. 



termination of some of the parameters including E$i, the 
thermal conduction efficiency, and the physics we include in 
our treatment of SN explosions. In addition our simulations 
do not contain a detailed treatment of metals. We demon- 
strate the effects of changing each of these parameters on 
the large scale properties of simulated quiescent disks. 

The value of E$i may differ greatly from unity, for ex- 
ample due to radiative cooling of the SN remnant. We in- 
vestigate the effect of moving Em away from unity and also 
look at changing the physics included in our analytic model 
for blast wave evolution, firstly by extending the simple Se- 
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Figure 19. Stellar density profile in the base rotating collapse 
simulation. Solid lines represent the best fit exponential and a 
de Vaucouleurs profile. The scale radii used in the exponential 
and de Vacouleurs fits are and rj, respectively. It is clear that 
even starting from an initial condition far from equilibrium we 
generate a stellar disk with a surface density profile similar to 
that in observed galaxies. 



Figure 17. Plot of the distribution of gas atomic in a galactic 
disk after 500Myr of evolution. Colours and plot dimensions are 
matched to those in Levine et al (2006) for easy comparison to 
observations. The inner circle represents the radius from the cen- 
tre of the galaxy to the position of the sun. The outer radius is 
where the observations of Levine et al (2006) are truncated. The 
simulated MW has a surface density profile in close agreement 
with that observed by Levine et al. 
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Figure 18. Density profiles for the ambient material (grey lines) 
and the molecular clouds (black lines). The solid lines represent 
the low resolution rotating collapse simulations and the dotted 
lines represent the highest resolution simulations. Agreement be- 
tween the high and low resolution simulations is very good. 

dov so lution with the fitting formula due to iTang fc Wanel 
l|2005l ) and secondly by investigating the effects of prevent- 
ing radiative cooling in supernova remnants. 

Altering the value of £51 has, as expected, two main ef- 
fects. Firstly an increased supernova efficiency can eject gas 
from the galactic disk more efficiently and quenches star for- 
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Figure 20. Fraction of gas in the molecular phase as a function 
of distance from the centre of the galactic disk for the rotating 
collapse model at two different times. The solid line represents the 
same data as plotted from M31 (Dame et al (1993)), with the x- 
axis representing position along the major axis of the galaxy. The 
simulated galactic disk is in good agreement with observation. 



mation very quickly . Secondly, the gas disk in simulations 
with higher supernova efficiency is found to be less centrally 
concentrated; the reverse is also true in simulations with a 
low supernova efficiency. This is a result of supernova feed- 
back ejecting gas from the galactic disk more efficiently in 
the runs with a high value of £51. We find that the fiducial 
value of £51 = 1.0 provides a good match with the observed 
properties of disk galaxies. 

More accurate modelli ng of the evolu t ion of supernova 
blast waves (by using the ITang fc Wangl (|2005l ) fit to the 
blast wave evolution) does not significantly change the prop- 
erties of the galaxy. Assuming that a typical supernova rem- 
nant expands for approximately ~0.3 Myr before being dis- 
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Figure 22. The relation between the pressure in the midplane of 
the galactic disk and the ratio of gas in the molecular phase to the 
atomic phase, R mo i- The dotted line represents the observational 
result (Eq J27))due to Blitz & Rosolowsky (2006). The simulated 
galaxy is in good agreement with observation. 



persed, that the mean ambient temperature of the ISM is 
10 6 K and the mean density is 10~ 2 atoms per cm 3 then 
the difference in the blast wave radius between the pure Se- 
dov solution (1 101) and the modified fitting formula due to 
iTang fc Wand (120051 1 (Eq JTTJ) is never more than 20%. 
This is demonstrated in Fig [2] which shows a pure Sedov 
blast wave compared with the simulated result from an SPH 
simulation of a blast wave in a hot (10 6 if) medium. 

The net result of having larger blast waves is that the 
porosity of the ISM increases at a greater rate and the delay 
between a supernova explosion occurring and its thermal 
energy being injected into the ambient medium is decreased. 

Finally we can switch off the radiative cooling in super- 
nova remnants. The effects here are much more dramatic. In 
the dense galactic disk a supernova remnant can typically 
radiate away over 90% of its thermal energy before being 
disrupted. Switching off radiative cooling in supernova rem- 
nants, therefore, has a very severe effect in our galaxy sim- 
ulations, and leads to the almost immediate suppression of 
star formation and much of the material in the galactic disk 
is ejected from the galaxy. This demonstrates that radiative 
losses from supernova remnants are of crucial importance in 
our models. 

Our simulation code does not contain a detailed descrip- 
tion of mass feedback from supernovae. Therefore we need 
to verify that the expected evolution in metallicity over the 
timescale of the simulation will not substantially affect the 
properties of t he galactic disk. 

Following lHarfst. Theis fc Henslerl |2006) we use a sim- 
ple analytic model to estimate the change in metallicity over 
the timescale of a typical simulation then run simulations 
with metallicites bracketing this range. By assuming that 
stars form at a constant rate of 1M@ /yr, that there is one 
type 2 supernova event per 125Mp) of stars form ed, and 
using metal yields due to IWooslev fc Weaver! (|l995l ) we es- 
timate that over the lGyr timespan of one of the quiescent 
disk simulations the average metallicity of the galaxy should 
not change by more than O.O4Z0. 

The base simulations have a metallicity of I.QZq, two 



additional simulations were run with metallicities of 0.5Zq 
and 1.5Zq, far outside the metallicity evolution range ex- 
pected in our quiescent disks. The total amount of stars af- 
ter lGyr in the high metallicity run is 5% higher than in the 
base run. The low metallicity run contains approximately 
5% less stars than the base run. This trend arises because 
the radiative cooling rate of the ambient phase is related to 
the ambient gas metallicity. All the properties of the three 
simulated disks agree to within 10% with the properties of 
the base simulations. The reason for the surprisingly small 
dependence of the galaxy properties on the metallicity of 
the gas is that in the quiescent disk the gas is maintained 
at a temperature of ~ 10 4 K. At this temperature the ra- 
diative cooling function does not change very much with 
metallicity. We do, however, note that a full treatment of 
the metal enrichment of the gas is necessary in perform- 
ing fully cosmological simulations since here the metallicity 
of the gas may significantly a f fect t he way in which it col- 
lapses (e.g. IScannapieco et all (|2005l 0. In these simulations 
we have additionally ne glected the change in stellar lifetim e 
with metallicity (see e.g. lRaiteri. Villata fc Navarrol (|l996t) ). 
Since we are simulating a quiescent disk, it is not expected 
that changes to the lifetimes of massive stars will have a 
large effect on the propert ies of the galactic disk. 

However, as noted by lHarfst. Theis fc Henslerl ([2006) a 
detailed prescription for the yields from supernovae is nec- 
essary if we want to simulate the early evolution of a galaxy. 

One of the most poorly constrained parameters is the ef- 
ficiency of evaporation of molecular clouds through thermal 
conduction. Magnetic fields and turbulence may affect the 
amount of thermal conduction by a large amount (M077). 
We ran quiescent disk simulations with the efficiency of ther- 
mal conduction moved by a factor of five in each direction 
(GAL_BASE_LOCON; GAL_BASE_HICON). More efficient 
thermal conduction leads to a lower density of molecular 
clouds in the galactic disk, as well as making the cloud mass 
function more shallow. The star formation rates are affected 
by a similar amount, in the simulations with a high ther- 
mal conduction rate the star formation rate is depressed by 
an order of magnitude. As discussed in section [5721 the base 
value for the thermal conduction efficiency reproduces many 
of the observed properties of the MW. 



6 CONCLUSIONS 

Motivated by the fact that we cannot reasonably resolve 
the Jeans scale for molecular clouds in galaxy simulations 
we have introduced a new star formation and feedback pre- 
scription. We model the ambient phase of the ISM using a 
hydrodynamic simulation code and the unresolved molecu- 
lar gas using a sticky particle prescription. Our model leads 
to a tightly self-regulating multiphase ISM. The multiphase 
nature of our star formation prescription avoids a lot of 
the problems of overcooling that were present in the first 
generation of star formation models. With the exception of 
the parameter that controls the molecular cloud coagula- 
tion timescale, v s tick, all the parameters in our model can 
be tightly constrained by observation, leaving the cloud co- 
agulation timescale as a free parameter that we can adjust 
to match the observed properties of galaxies. Where pos- 
sible our model of the ISM has been formulated in such a 
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way that the results of a simulation should be independent of 
mass resolution. We demonstrated that the large scale prop- 
erties of our simulations were unchanged over a two orders 
of magnitude shift in mass resolution. 

We have applied the sticky particle star formation 
model to three different types of simulation: a simple one 
zone model, the rotating collapse of a gas and dark matter 
sphere and a model of a quiescent galactic disk. After using 
the one zone simulation to set the value of the parameters 
that cannot be determined observationally, the sticky par- 
ticle model can be applied to the other simulations without 
any parameter changes. 

The simulations of a quiescent disk galaxy reproduce 
the observed Schmidt law with a slope of 1.4 due to the 
opposing effects of cloud coagulation and feedback effects. 
The galaxy also developed a natural three component ISM. 
Finally we observe supernova heated gas in the galaxy being 
ejected from the disk either in the form of a galactic fountain, 
or, when the star formation rate (and associated supernova 
rate) is sufficient, in the form of strong bipolar outflows. 
Both of these results arise as a natural consequence of the 
physics included in our star formation prescription. 

Simulations of the collapse of a rotating sphere of dark 
matter and gas reproduced, many of the observed proper- 
ties of galactic disks, beginning from an initial condition 
well out of equilibrium. In particular we observe a stellar 
disk with distinct bulge and disk components, well fitted by 
the standard exponential and de Vacouleurs density profiles. 
The fraction of molecular gas in the disk as a function of ra- 
dius is reproduced, and agrees well with recent observations 
of nearby galaxies. The observed relation between the disk 
midplane pressure and the fraction of molecular clouds is 
also reproduced. We also observe star formation rates com- 
parable to those in disk galaxies and note that our model 
reproduces the formation of stars in the spiral arms of the 
galaxy. 

Our preferred values for most of the parameters are dis- 
cussed in section 2] and are usually constrained by observa- 
tion. We varied the values of those parameters that are only 
weakly constrained and foind that our preferred value of 
£51 = 1 find that we obtain reasonable results. Including 
different physical descriptions of the evolution of supernova 
blast waves makes only a modest difference to our results. 
Finally we note that in the quiescent disk simulation our re- 
sults depend relatively weakly on metallicity; although this 
will not be the case in fully cosmological simulations, for 
which a full treatment of metal enrichment via supernova 
feedback is necessary. 

A natural continuation of this work is to extend our in- 
vestigations to higher redshift through the use of fully cos- 
mological simulations, and to explore the behaviour of the 
ISM in colliding galaxies. This work is currently being pur- 
sued. 
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APPENDIX A: LIST OF SYMBOLS 

a c : Slope of the molecular cloud mass-radius relation. Eq[2] 

Ch : Sound speed of the ambient gas phase 

e* : Fraction of a GMC converted to stars in a collapse. 

r\ : Fraction of cloud velocity lost to 'cooling' collision 

i?5i : Energy ejected per SnII in units of 10 51 ergs 

Eb ■ Total energy in a supernova blast wave 

Em : Total kinetic energy in molecular clouds of mass m in 

a given volume 

/d : Filling factor of cold clouds 

/m(o"i,a2) : Fraction of collisions between clouds with 

velocity dispersions oi and ui that lead to mergers 

K(m, m') : The kernel for aggregation of clouds of masses 
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m and m' . Eq. [3] 

A : Constant of proportionality relating cloud mass and 

destruction rate by thermal conduction. Eq. [22] 

An : Normalised radiative cooling rate 

Anet : Net radiative cooling rate (ergs cm~ 3 s _1 ) 

M c : Mass of a molecular cloud 

Af c h a r : Characteristic mass resolution of a simulation 
M re f : Reference cold cloud mass. Eq[2] 
Mj^mjn : Minimum allowed star mass 
Af* jmax : Maximum allowed star mass 

rib : Density internal to a supernova remnant in atoms / 
cm 3 n c : Density of a molecular cloud in atoms / cm 3 
rih : Density of the ambient medium in atoms / cm 3 
A^sf : The slope of the schmidt law. Eq[8] 
n(m, t) : The number of clouds with masses between m and 
m + dm 

N(m, t) : The number density of clouds with masses 
between m and m + dm 

cf> : Efficiency of destruction of cold clouds by thermal 
conduction 

Q : Porosity of the interstellar medium. Sec. 13.4.21 

r c : Radius of a molecular cloud 

r rc f : Reference cold cloud radius. Eq[2] 

rb : The radius of a spherical blast wave 

p c : Mean density of molecular clouds contained in a volume 

ph : Mean density of ambient gas contained in a volume 

pth '■ Density at which ambient gas becomes thermally 

unstable 

Psfr '■ Volume density of star formation 

Tb : Mean temperature inside of a supernova remnant 

T c : Internal temperature of cold clouds 

Th : Temperature of the ambient gas phase 

itb : Thermal energy per unit mass of supernova remnants 

7i c : Thermal energy per unit mass of the cold clouds 

ith : Thermal energy per unit mass of the ambient phase 

E : Cross section for collision between clouds. Eq. U 

Econd : Efficiency of thermal conduction. Eq 1181 

«a PP : Relative approach velocity of two molecular clouds 

■"stick ' Maximum relative velocity for cloud merger 

x : Slope of the stellar IMF 



and rri2 that lead to mergers is given by 



/m(oi, 02) = 



I Vl + ^tick N r/ V l — ^Stick 

erf ( — ] — erf 



(B2) 



Using this definition of / m the Smoluchowski equation 
becomes 



dn _J_ 

~dt ~ZV 



n{m , t)n(m — m ,t)K{m , m — m ) 



n(m, t) 



n(m ,t)K(m,m')fm{o- m ,a m i)dm . 



(B3) 



As discussed in section [3^21 clouds of mass m may gain 
or lose kinetic energy in three ways: clouds of mass m' and 
m — m' may merge to form extra clouds of mass m. Clouds 
of mass m may merge with clouds of any other mass to 
decrease the number of clouds of mass m. Finally clouds of 
mass m can interact gravitationally with any other clouds, 
thus losing kinetic energy. These three processes are termed 
gain, loss and cooling. 

Gain processes may be represented in the following way, 
where we have integrated over m' such that the two particles 
that merge have masses that sum to m 



-E'gain 



00 rvi — 00 rV2 — vi + i> s tick 



J v-i — — 00 J — V 



2=«1 -"stick 



P(vi)P(v 2 )n(m',t) 



n(m — m', t)K(m', m — m')f m (m', m — m')£J/j 
dv2dvidm . (B4) 



APPENDIX B: ENERGY TRANSFER 
THROUGH COAGULATION 

Our formalism to treat the evolution of the mass function of 
clouds internal to each of the 'multiple cloud' particles will 
start from the Smoluc howski equation of kinetic aggregation 
l|Smoluchowskil (119161 )1 



dn 

~M ~ 2V 



— I n(m',t)n(m — m',t)K(m' ,m — m)dm 
* Jo 



V 



n(m,t) I n(m' ,t)K(m,m')dm , (Bl) 



where n(m, t) represents the number of clouds with masses 
between m and m + dm contained within a volume, V and 
K(m, m') represents the kernel for aggregation of clouds 
with masses m and m', as defined by Eq @. 

The fraction of collisions between clouds of masses mi 



P(vi) and P{vi) are the probability distributions velocities 
Ml and V2 and are assumed to be gaussian with standard 
deviation <ti and 02 respectively. , Ef represents the final 
kinetic energy of a collision between particles of masses m' 
and m — m' . Ef is evaluated by considering conservation of 
momentum, 



E f 



1 (m'vi + (m — m')v2) 2 



(B5) 



Eq (|B4|I then becomes 



Ee. 



n(m', t) 



f j\ roci 



/ n( 
Jo 

/oo rvi+v stick 2 
-co J vi — " st i c k ^m'&m — m 

1 / (m'vi + (m — m')v2) 2 



m — m ,i)K{m , m — m) 



dv2dv\dm . 



(B6) 
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and Eq ()B6|I may be written 



Sgain - ^ 



n(m',t) / x 



n(m — m ,t)K(m , m — m ) 



oo ^ m ^ m — m 



r i+ "" a e -( ^_ m J l/ (m' 8l + (r»-m> 2 ) 2 

•'"l -"stick 



dvidv\dm! 



(B7) 



APPENDIX C: THE SOLUTION OF THE 
COAGULATION EQUATIONS 

In our simulations we solve the discrete versions of Eq (|B3[) , 
Eq ((5]l and Eq ©, where by assuming that cloud mass is 
quantised into N bins characterised by an index, i, where 
Mi = iMo we can write 



1 71 j"*«*->i 

= 2V 7 K nfij n i n i -yYl K ikf^nj, (CI) 

i+j — fc j — 1 



The total kinetic energy of particles of mass m may also 
be decreased by mergers between particles of mass m and 
any other mass (the second process in the list). Similarly to 
Eq (|B6[) . the rate of energy loss may be written 

• n{m,t) 



2tt 



1 



r oo 

/ n(m\ t)K(m, m') 

JO J - oc ^m^m' 



«l+"stick 



mv'f 



dv2dvidm' . 



(B8) 



Finally, the total energy of particles with mass m may be de- 
creased by collisions between particles of mass m and parti- 
cles of any other mass that occur at relative velocities greater 
than u s tick- In this case, the velocity of both particles is de- 
creased by a factor 77 relative to the centre of mass. For a 
collision between particles of masses mi and m 2 (velocities 
vi and V2) the final velocity of particle 1 (denoted v[) is 
evaluated by conservation of momentum 

Vl =77(l>l - Hcom) + fcoin 

mivi + m-zvi 



A TTi 1 12 1 2 
AE =-rriiv 1 — —m\V 1 . 



(B9) 



Using these definitions, the change in energy of a particle 
of mass mi by gravitational cooling with a particle of mass 
777,2, denoted e is given by 



mi \ 2 



\yx{\ - a 2 ) - vj(3 2 - viv 2 a/3^ . 



where 



a = 77 + 



m 1 



mi + m,2 
m2 



Till + 7712 



(BIO) 

(Bll) 
(B12) 



In a similar way to Eq (|B6|) the energy loss via this process 
may be written 

■ n(m,t) f°° . , . 
E coo i = — - / n(m ,t)K(m,m ) 



2tt 



/ 



"1 — "2 1 >"stick 



" ( V2v m , ) 



dv2dvidm' . (B13) 



E k — -S'gain E\ oss E coo \ , 



Ck = 



E k - \a\M k h k 



(C2) 



(C3) 



o k n k M k 

where subscripts represent different mass bins. Kij = 
K(Mi, Mj) = K(iMo, kMo). The superscript m represents 
that / is a cross section for particle mergers 

To demonstrate the technique for solving these equa- 
tions we will consider the numerical solution of the simple 
Smoluchowski (Eq ([BlJ), which when written in a discrete 
form takes on the following form 

. 2-1 AT 

hi = n 3 n i-3 K ij ~ y njKij. (C4) 

3=1 3=1 

Following iBenson. Kamionkowski fc Hassan! (|2005l ) Eq 
(|C4|) can be rewritten in the form of a matrix equation 

n = B k, (C5) 

the vector k has N x N elements corresponding to 
K(mi,mj). The kernel matrix, B has N x N x N elements 
and may be written more explicitly as 



Bij k kj k , 

jk 



where 



B 



ijk 



rijn k 
V 



(C6) 



(C7) 



S represents a Kronecker delta function. We solve Eq (|C5|l 
implicitly using an iterative method. 

The solution of the equations that govern energy ex- 
change between clouds (Eq lfB7)l. Eq JB8) and Eq (|B13[) ) 
is the same as for the solution of the Smoluchowski equation 
in that we will write the equations in the form of the linear 
multiplication of two matrices and then solve this equation 
implicitly. In order to simplify the notation in this section 
we will denote the terms in the three equations that are in- 
side of the integrals over velocity as £. Explicitly for the case 
of the equation for energy gain (Eq (|B7[) ): 



dv2dv\ 



C G (m,m') = — [°° e (^O 

27! " J -co a™<J m i 



L 



-e V 



(mvi + m'viY 



(C8) 



The corresponding terms in the equations for energy 
loss (|B8|) and cooling (|B13|) are denoted £ L (mi, 777,2) and 
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£ c (mi, 7712) respectively. Note that the definitions of £ in- 
clude the factors of 27r and - from throughout the equations. 

The equation for the total evolution of the energy of a 
system of coagulating and cooling particles may be written 
in terms of these new functions as: 

f-OO 

E(m) = / n(m',t)n(m — m',t)K(m',m — m')£ G (rri ,m — m')dm' 
Jo 

f-OO 

— n(m,t) / n(m' ,t)K(m,m')^ L (m,m')dm' 

Jo 

f-OO 

— n(m,t) / n(m' ,t)K(m,m')£ c (m,m')dm' (C9) 

Jo 

Which when discretized and rearranged becomes 

i-l N 

Ei = J2 runjK^. - n< ( ^ u,K %] + C )) (CIO) 
j=i j=i 

The subscripts represent different mass bins (jij = n(jMo)). 
Our goal is to rewrite Eq ()C10|) in the form of a linear 
multiplication of two matrices 

E = Ck, (Cll) 

where k is defined in the same way in the solution of the 
Smoluchowski equation, that is: kij = K{irii,mj). The form 
of Cijk that is consistent with Eq (|C 10|) is given by: 

Cijk = rijUu (br.j+kifk - Sik (£fk + C^fe)) (C12) 

This form for Cijk is functionally equivalent to (Eq 
(|C7I) ) so the solution may proceed in exactly the same way 
as for the Smoluchowski equation, the only difference is the 
form of the matrix B 

The calculation of the quantities £ G , £ c and £ L is com- 
putationally very expensive so they are initialised once into 
a lookup table at the start of every simulation and obtained 
by bilinear interpolation thereafter. 

This paper has been typeset from a TgX/ ETgX file prepared 
by the author. 



